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FOREWORD 

MAPSEP (Mission Analysis Program for Solar Electric Propulsion) 
is a computer program developed by Martin Marietta Aerospace, Denver 
Division, for the NASA Marshall Space Flight Center under Contract 
NAS8-29666. MAPSEP contains the basic modes: TOPSEP (trajectory 

generation), GODSEP (linear error analysis) and SIMSEP (simulation). 
These modes and their various options give the user sufficient flexi- 
bility to analyze any low thrust mission with respect to trajectory 
performance, guidance and navigation, and to provide meaningful sys- 
tem related requirements for the purpose of vehicle design. 

This volume is the first of three and contains the analytical or 
functional description of MAPSEP. Subsequent volumes relate to opera- 
tional usage and to program logical flow. 
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1. INTRODUCTION 

A major requirement for spacecraft systems design is the effec- 
tive analysis of performance errors and their impact on mission success. 
This requirement is especially necessary for low thrust* missions where 
thrust errors dominate all spacecraft error sources. Fast, accurate para- 
metric error analyses can only be performed by a computer program which 
is efficiently constructed, easy to use, flexible, and contains model- 
ing of all pertinent spacecraft and environmental processes. MAPSEP 
(Mission Analysis Program for Solar Electric Propulsion) is designed 
to meet these characteristics. It is intended to provide rapid evalu- 
ation of guidance, navigation and performance requirements to the 
degree necessary for spacecraft and mission design. 

The baseline design of MAPSEP was taken from a previous study 
effort (Reference 1). Suitable modifications to the design were made 
which reflected subsequent operational experience (References 2 and 3) 
and actual construction and testing of MAPSEP. Considerable knowledge 
was also gained in the construction and usage of the engineering ver- 
sion of MAPSEP which was actually three separate programs that corre- 
sponded to the modes (Figure 1-1): TOPSEP, GODSEP and SIMSEP. Driving 
considerations in the program design and construction were: realistic 

vehicle and environment modeling consistent with preliminary vehicle 
design, flexibility in usage, computational speed and accuracy, mini- 
mum core utilization (for turnaround time and operating costs) and 
maximum growth potential through modularity. 
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This document is the first of three volumes. Contained herein 
is a brief introduction to MAPSEP organization and detailed analyti- 
cal descriptions of all models and algorithms. These include, for 
example, trajectory and error covariance propagation methods, orbit 
determination processes, thrust modeling and trajectory correction 
(guidance) schemes. This analytic background is necessary to fully 
understand program operation and to maximize program capability with 
respect to the user* 

The second volume is a description of program usage, that is, 
input, output, recommended operating procedures and sample cases. 

The third volume is a detailed description of internal MAPSEP struc- 
ture including macrologic, variable definition, subroutines and logi 
cal flow. 

The Analytic, User's and Program Manuals represent a version 
distinct from the baseline interplanetary program. The baseline 
MAPSEP, delivered to NASA/MSFC in October, 1974, emphasized inter- 
planetary low thrust missions with limited Earth orbital capability. 

The current MAPSEP version, as reflected by the following documentation, 
concentrates specifically on the use of low thrust primary propulsion 
for Earth orbital applications. 

Many features of the interplanetary version were removed or made 
inoperable because of Earth orbital requirements. Some computational 
inefficiency still exists because of the limited contract scope. In 
any case, Earth orbital MAPSEP contains the capability of analyzing 
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almost any currently projected low thrust mission from low Earth 
orbit to "super" synchronous altitudes. Furthermore, MAPSEP is 
sufficiently flexible to incorporate extended dynamic models, 
alternate mission strategies, and almost any other system require- 
ment imposed by the user. 

As in the interplanetary version, Earth orbital MAPSEP 
represents a trade-off between precision modeling and computational 
speed consistent with defining necessary system requirements. It 
can be used in feasibility studies as well as in flight operational 
support. Pertinent operational constraints are available both 
implicitly and explicitly. However, the reader should be warned 
that because of program complexity, MAPSEP is only as good as the 
user and will quickly succumb to faulty user inputs, sometimes 
in an obvious fashion, and sometimes in a very subtle manner. 
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2. PROGRAM DESCRIPTION 

This section summarizes MAPSEP' s function and use, and MAPSEP 
structure. These areas are discussed in greater detail in the user's 
manual and programmer's manual, respectively. 

As mentioned earlier, MAPSEP is composed of three primary modes. 
Each mode is intended to serve a given function in the mission design 
sequence. TOPSEP (Targeting and Optimization for SEP) is used to 
generate numerically integrated trajectories consistent with dynamic 
and system constraints. Performance data and related sensitivities 
are computed in the process of trajectory generation but can also be 
obtained by parametric application of TOPSEP. Indeed, each mode 
readily lends itself to parametric use which is a necessary feature 
for mission and system design studies, 

GODSEP (Guidance and Orbit Determination for SEP) is used to 
perform a linear covariance analysis about a selected reference tra- 
jectory, generated by TOPSEP. Various dynamic and measurement related 
error sources are applied in a statistical sense to compute trajectory 
error covariances. These covariances, corresponding to estimation un- 
certainty (knowledge) and to actual trajectory deviations from the 
nominal (control), are propagated through a sequence of mission events: 
thruster switching, navigation measurement /state update, trajectory 
correction (guidance), etc. Thus, GODSEP computes a time history of 
the ensemble of all expected trajectory errors, and in the process 
displays such useful system parameters as required thrust control 
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authority, predicted terminal miss, additional fuel expenditures for 
off-nominal performance, etc. 

SIMSEP (Trajectory SIMulation of SEP) is used in the latter 
stages of system design. It deterministically simulates the trajec- 
tory including the application of discrete dynamic errors. Trajectory 
corrections are simulated in an operational sense through a thrust 
update design and execution process. Navigation is simulated by 
sampling estimation error covariances (generated by GODSEP) prior to 
each guidance event. By operating SIMSEP in a Monte Carlo fashion, 
any desired number of simulated missions can be obtained, and estimated 
singly or collectively in statistical displays. 

A fourth "mode", REFSEP (REFerence SEP), is actually an expansion 
of TOPSEP to provide a greater amount of trajectory and navigation 
related data for a particular reference trajectory. 

Each. mode of MAPSEP uses a common trajectory propagation routine, 
TRAJ. This guarantees trajectory reproducibility among modes, TRAJ 
integrates the equations of motion in Encke form using a fourth order 
Runge-Kutta scheme. Covariance propagation and transition matrices 
are computed by integrating variational equations simultaneously with 
the equations of motion. An option exists in GODSEP which stores a 
complete set of trajectory parameters and transition matrices as it 
is generated by TRAJ onto a magnetic disc called the STM file so that 
subsequent error analyses will not have to regenerate the data. The 
information on disc can then be transferred to magnetic tape for per- 
manent storage, if desired. A more limited option is available for 
TOPSEP and SIMSEP which stores only the initial (input) trajectory 


data on disc. 
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Input to MAPSEP is primarily through cards using the NAMELIST 
feature, with supplementary means depending upon mode and function 
(Table 2-1). All modes require the $TRAJ namelist which defines 


Mode 

INPUT 

OUTPUT 

Namelist Formated Tape 

Cards (or disc) 

Punched Tape 

Cards (or disc) 

TOPSEP 

$TRAJ None STM 

$T0PSEP 

None STM 

GAIN 

GODSEP 

$TRAJ 

$G0DSEP Event STM 

$GEVENT Data GAIN 

States 

Covariances STM 

Guidance GAIN 

SUMARY 

SIMSEP 

JTRAJ None STM 

ISIMSEP 

fGUID 

Statistics STM 

GAIN 
SUMARY 

REFSEP 

$TRAJ Print STM 

Events 



TABLE 2-1. MAPSEP Input/Output 

the nominal trajectory and subsequent mode usage. However, if 
recycling or case stacking is performed it is not necessary to in- 
put $TRAJ again unless desired. The second namelist required for 
each mode corresponds to mode peculiar input and bears the name of 
that particular mode. Additional namelist, formated cards, and tape 
input are generally optional. Besides the standard printout asso- 
ciated with MAPSEP, auxiliary output can be obtained which will 


facilitate subsequent runs. 
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The structure of MAPSEP is organized into three levels of 
"overlays" which are designed to minimize total computer storage. 

At any given time, only those routines which are in active use are 
loaded into the working core of the computer. The main overlay 
(Figure 2-1) is always in core and contains the main executive, 
MAPSEP, and all utility routines that are common to the three modes. 
The primary overlays contain key operating routines of each mode, 
that is, those routines which are always needed when that particular 
mode is in use. Also included as a primary overlay is the data 
initialization routine, DATAM, where $TRAJ namelist is read, trajec- 
tory and preliminary mode parameters are initialized, and appropriate 
parameters are printed out. 

The secondary overlays contain routines which perform various 
computations during a particular operational sequence. Included are 
data initialization routines, analgous to DATAM, which operate on 
mode peculiar input and perform mode initialization. An example of 
core usage in the changing overlay structure may be provided by a 
standard error, analysis event sequence. Error analysis initializa- 
tion is performed by the overlay DATAG. Transition matrices are then 
read from the STM file, the state covariance is propagated to a 
measurement event, and the overlay MEAS is called, which physically 
replaces, or overlays, the same core used previously by DATAG. 
Similarly at a guidance event, overlay TRAJ will replace MEAS to 
compute target sensitivity matrices and overlay GUID will then 



MAIN 


MAPSEP 


PRIMARY 


SECONDARY 



Figure 2-1. OVER 




•~ | utility"] 



STRUCTURE 
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replace TRAJ to compute guidance corrections. Overlay switching is 
performed internally and is transparent to the user. 

All of the routines and structure of MAPSEP are constructed to 
minimize core storage (thus reducing turn-around time and computer 
run cost) yet retain the flexibility needed for broad analysis re- 
quirements. Furthermore, routines are built as modular as possible 
to reduce the difficulties in future modifications and extensions. 



3 . NOMENCLATURE 


The following symbols are used throughout the Analytic Manual, and 
to a great extent in the User f s and Program Manuals, However, deviations 
from these symbols may occur in localized discussion if required for pur- 
poses of clarity* 


SYMBOL 


a 

c 



E 

F 

g 


H 

I 

K 

m 

P 


P 

o 

Q 

r 

s 


S 


DEFINITION 

propulsive acceleration 

propulsive exhaust velocity 

cross covariance between i and j parameters 

target error index 

dynamic variation matrix 

performance gradient or thrust transfor- 
mation matrix 

observation sensitivity matrix (WRT state) 
identity matrix 
filter gain matrix 
spacecraft mass 

covariance of state deviations or elec- 
trical power or projection operator 

effective power at 1 AU 

dynamic (thrust) noise matrix 

spacecraft position 

solve-for parameters 

target sensitivity matrix WRT control para- 
meters 




U 

t 

SYMBOL 

DEFINITION 


T 

event time, or target variables, or 
thrust 

. 

u 

dynamic (consider) parameters 


ut 

thrust acceleration proportionality (throttling) 


U 

V 

Control parameters 

spacecraft velocity or measurement para- 



meters 


w 

ignore parameters 


W 

weighting matrix 


X 

spacecraft state 


r 

guidance matrix 


1 

propulsive efficiency 


e 

transition matrix of dynamic parameters 



gravitational constant 



relative range 


a 

standard deviation 


r 

correlation time of thrust noise 



transition matrix of augmented state 



state transition matrix from time t^ to 
fc k+l 


60 

thrust noise 


SUBSCRIPT 

DEFINITION 



state control covariance 


( 'k+l,k 

matrix evaluated over time interval t^ 

to s+i 


< >k 

state knowledge covariance 



12 


SUBSCRIPT 


DEFINITION 


< >v 

< >. 
< >* 


planet related parameters 
solve-for parameters 
measurement consider parameters 
ignore parameters 
spacecraft -state parameters 


MISCELLANEOUS DEFINITION 


G&N 

guidance and navigation 

OD 

orbit determination 

PGM 

Projected Gradient Method 

S/C 

spacecraft 

SEP 

solar electric propulsion 

WRT 

with respect to 

E[ ] 

expected value operation 

( ) + 

post-event value 

( )“ 

pre-event value 

• 

( ) 

time derivative 

A 

( ) 

unit vector 
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4.0 TRAJ 

Essential to any program used for performance and navigation 
analysis is an accurate, but computationally efficient, trajectory 
propagation routine. This routine must contain realistic models of 
the dynamic processes acting on and performed by the spacecraft. In 
MAPSEP, the subroutine TRAJ fulfills this role. TRAJ is designed to be 
used by the three modes TOPSEP, GODSEP and SIMSEP, and is capable 
of duplicating the same trajectory in all modes. 

The trajectory overlay TRAJ propagates low thrust, or ballistic, 
trajectories relative to Earth, using Encke's formulation of 
the equations of motion, from any epoch to a termination condition. 

TRAJ can optionally propagate the state covariance or the state 
transition matrix along the trajectory for either the basic state 
or an augmented state. Two of the most important features incorpo- 
rated into TRAJ are the variable integration step algorithm and 

trajectory repea tibi li ty . 

4,1 Equations of Motion 

In Encke's formulation, all accelerations other than those due 
to the gravity of the Earth are called perturbing accelerations. 

a 

Position (r ) and velocity (r ) vectors are computed relative of this 
"~C c 

primary body using two body formula. The deviation vectors from 
the reference conic position and velocity vectors, £ r_ and & r_, 
respectively, are the direct results of numerically integrating the 
sum of the perturbing accelerations. The true position and velocity 
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vectors are, respectively, 

r = r + Sr. 

— ”C 

4 * • 

r = r + £ £ 


Since r and r can easily be computed from conic formula (See 
— c — c 

Appendix 1), the only problem is to compute r and r. 

Let £r be the acceleration deviation from two body motion, such 
that Sr. is the sum of a11 Perturbing accelerations, e.g., other 
bodies, thrust, etc. 

N 

i + Si] - [r +f(« t ) r.] 

c i=l ^ 

+ l + %+i J2 


The first term is the difference between two body and perturbed two 
body motion. The second term is the sum of the perturbing accelera- 
tions due to the sun and/or moon. The third term (a) is the accelera 
tion due to thrust. The fourth term (a^) is the acceleration due to 
radiation pressure, and the fifth term, the acceleration resulting 
from a nonspherical mass distribution in the Earth. 

The first term is computed from the following equations (See 
Reference 4), 


r , *<3 + 3* + */) 

£(°0 = T /2 

1 + (1 + cO 


v = ( Sn - 2r) * Sr 

2 

r 
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where r. is the true s/c position vector relative to the primary body 
and ja is the gravitational constant of the primary body. 

To compute the second term, the heliocentric position vectors, 
r., of the perturbing bodies, are computed from mean analytical 
orbital elements to obtain 


^ = r + r p - r. 


f(rf-) = 


*i = 


^i 



3 +3 ©t ± +°( 
1 + (1 +*.) 



- 



2 r 


9i 


r «i 




where 

The 


is Che heliocentric position vector of the Earth 
acceleration vector due to radiation pressure is 


% 


1.024 x 10 8 A C r 


2 

tn r 


A 

r 


where 1*024 X 10^ 
m 
r 
A 
C 

r 

The option exists 


- solar flux constant 

- instantaneous mass of the s/c 

- heliocentric position vector of the s/c 

- effective cross sectional area of the s/c 

- coeffient of reflectivity. 

in TRAJ, to include or exclude the effects of 
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\ 


radiation pressure when propagating both planetary and interplanetary 
trajectories* Before defining the acceleration due to thrust ji, we 
will model the power subsystem* There are two power subsystem models 
used by TRAJ for low thrust* They are solar electric and nuclear electric 
The power to the thrusters (p) is 


P G , 
o o 




exp 




solar 

electric 


P = 


P if P > P or r < r . 
max, max mm 


solar 

electric 


P ° “P [- P l (t - t UL> ] - V 


nuclear 

electric 


P - Power available (at 1 AU for solar, at energization for 
u nuclear electric 

- (Empirical) Constants defining solar array characteristics 

r - Heliocentric position magnitude of the S/C 

P^ - Power decay constant 

t - Time from epoch 

t - Time delay 

JJi-i 

P - Housekeeping power 

HK 

P - Maximum allowable solar electric power 

max 


r . - Heliocentric distance at which P reaches P 

min max 

The exponential term in the solar electric expression describes the 

degradation of the solar array as a function of time. 
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C is a scale factor which represents degradation of the solar 
o 

cells due to particular bombardment in the Earth's radiation zones. 
The radiation field is assumed to be axially symmetric about the 
magnetic pole. Power degradation effects are modeled after 
Reference 9 where the power loss factor is defined by 

-= o (lo 8lo F) 10 

C ~ e 
o 

where a is a constant and F is the cumulative particular flux 
o 

(measured in equivalent 1 MEV electrons per square cm)- F is eval- 
uated by integration of the flux rate, 


F * a io e 


V 7 a. £■ 

JLi 1 1 

2 =1 


where a^, a^, ...» are empirically derived constants and the f^' s 

are functions of vehicle distance, R, in Earth radii, and geomagnetic 
latitude, S . This flux model is most effective above a 6000 km altitude 

cos S 


f l = 


cos 


R 

cos £ 


cos £ sin £ 


R 


£ 4 " 


f 5 = 


f 6 = 


sin & 


sin^S 


,3 

sin S 


R 

. 3 


R 


f 7 = 


f 8 = 


f 9 


R 1/2 

COS S 

R 1 ^ 
cos 5 

. 1/5 


For numerical purposes, MAPSEP input and internal TRAJ computations 

14 

assume the flux units to be in 10 MEV electron per square cm. 

The thrusters provide the spacecraft with the ability to maneu- 
ver. By changing the orientation of the thrus t vector and maintain- 
ing that orientation for a long enough time, it is possible to steer 
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the spacecraft and "shape" the trajectory. The thrust controls are 
defined in terms of constant parameters over a given time segment of 
the trajectory. Thus the user or mode can specify the following 
controls for each segment (up to 40) : -* 

1. Thrust policy: Inertial pitch and yaw angles, In and 

Out of Plane Angles (orbit plane coordinates), or coast. 

2. Segment end time (referenced to launch) of the current 
segment. 

3. Throttling level, u^,. 

4. Pitch angle (or In Plane Angle), a^. * 

5. Coefficient, a^. 

6. Coefficient, a 2 * See Page 17 

V for definition 

7. Yaw Angle (or Out of Plane Angle), a^ A of a Q , ..., a $ 

8. Number of operating thrusters. 

9. Coefficient, a 

4 

10. Coefficient, a^ J 

The thrust policy is merely thrusting or not thrusting (coasting). 

Imposed coast periods result when the S/C passes into the shadow of 
the Earth. A complete discussion of the shadow model may be found in 
Appendix 9.8. During thrust, the user has an option as to the reference 
system for the acceleration vector. The two reference systems are 
inertial pitch and yaw Angles, Figure (4-1), and In and Out of Plane 
Angles, Figure (4-2). 

In addition, there are two types of pitch and yaw control policies 
which can be used as a function of the orbital mission type and desired 
controllability. Thus, the three available thrust policies are 
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0 = pitch = a + a t 4- a sin (E + 0- ) 

P F 0 1 2 1 


0 = yaw 

y 


= t + sin (E + 0 ^) 


OR 


I inertial 
pitch/yaw 


0 = pitch - a n + a 1 sin (E + a n ) 

p 0 12 

0 = yaw - a 0 + a. sin (E + a r ) 

y 3 4 5 


OR 


Y - in plane = t + sin 

£ = out of plane = a^ + t + 


(E + 0 t ) 
sin (E + # 2 ) 


I orbit 
I plane 


where t is the time from the start of the current control phase, 

E is the orbital eccentric anomaly, and and are phasing angles. 
For near circular orbits, use of the eccentric anomaly can cause 
difficulty because of the rapid motion of periapsis* MAPSEP thus 
has the option of using an alternate form of the fast variable E # 
Instead of the eccentric anomaly, a time dependent anomaly is used 

which is measured from the start of the current thrust phase. 
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PITCH 



X-Z plane 
is ecliptic 
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The acceleration due to thrust is 


I a | = 


0.002 ^ P u T 

m c 


.002 - Conversion factor 

rj - Averaged efficiency of hhe thrusters 

P - Power delivered to the thrusters 

u - Throttling level 
T 


m - Spacecraft mass 
c - Exhaust velocity 

The mass is numerically integrated from the equation 


m 


- ma 
c 


The thrust acceleration vector a/relative to the spacecraft for 
the two reference systems is 


a' 

= a 

cos 

e 

cos G 

X 



y 

P 

a * 

= a 

sin 

e 


y 



y 


a' 

= -a 

cos 

0 

sin 0 


z y P 

for the pitch-yaw system and 


a 1 

X 

- a 

cos S 

cos 

V 

a' 

y 

- a 

cos $ 

sin 

T 


a sin % 
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for the In and Out of Plane system. In order to transform a/ into a, 
which is in heliocentric ecliptic coordinates, we must define the 
matrix A such that 


where 


a = A a 1 


A = £x' ; Y ' ; Z ' j 

and Z' = -r/fr | 



— is the heliocentric spacecraft position vector and Z is the North 

s 

ecliptic unit vector for the inertial pitch/yaw System. For the 
In/Out of Plane system, A is defined in terms of the position and 
velocity vectors relative to the primary body. Let r be the 
position vector and v be the velocity vector, then A can be 
defined as 


A -[Aj 



where 


-1 " V / I V I 


r x v 
Ir x v| 


*2 = 


A 

“3 


x A^ 


and 
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Since the trajectory is made up of segments, TRAJ propagates 
the trajectory using the appropriate set of controls over the tra- 
jectory interval where they are in effect. Updates are automatically 
performed at the beginning of each new segment. 

The last perturbing acceleration, which is yet to be computed, 
is the acceleration resulting from a nonspherical mass distribution 
in the Earth. For TRAJ, several simplifying assumptions have been 
imposed such that only the effect of the J 2 harmonic coefficient 
included in the dynamic model. To begin, it is assumed that the 
Earth is axially symmetric. With this assumption, the gravitational 
potential function can be written as 

[ CO K 

1 - 2 J k ( P k (sin0) 

where V(r, 0 ) is written in terms of spherical coordinates r and 0 . 

Physically, r is the magnitude of the geocentric radius vector, and 

9 corresponds to the geocentric latitude measured relative to the 

equatorial plane. The coefficients (k = 2, 3, 4, . . . .) are 

harmonic coefficients which have been empirically evaluated for the 

Earth by satellite observations. R 0 q is the equatorial radius, and 

P (sin © ) is the k th Legendre polynomial of the first kind. The 
k 

expression for V(r,0) implicitly assumes that an equatorial 
inertial coordinate system is being used* In TRAJ where all dynamical 
calculations are performed relative to the ecliptic reference system, 
it is necessary to make appropriate coordinate transformations both 
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before the evaluation of V (r, 0), and after the desired acceleration 
has been computed. 

Expanding V(r,9) and neglecting all third and higher order terms, 
the gravitational potential becomes 

V(r,0 ) = 7 " r J 2 ( "^ ) P 2 (sin 

where (sin©) = k (3 sin I 2 0 - 1). The gravitational acceleration 
due to the nonsphericity is computed as a gradient of the potential 
function, i.e.. 


a j 2 = V V (r , 0 ) 


?v a , l 3v 

T7 r+ 7?i 


A 

a 


A 

+ (o) <J> 


where r, 0 and $ are orthronormal vectors defined in a spherical 
coordinate system. Performing the indicated differentiations and 
simplifying yields the acceleration acting on a S/C due to the 
attraction of an oblate planet as 



f -/* 3 f* J 2 R eq 

= L 7 


P^ (sin 0) 




3 JU J R sin© cos 9 
• 2 _e£ 


] 


0 


I 

The vector must be transformed from its spherical coordinate 
representation into a Cartesian frame. This is accomplished by the^ 

operation 


I 



where A is given by 
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x/r 

-y / > 

/ 2 2* 
/ x + y 

-xz/r Jx 2 + y 2 

A = 

y/r 

+x/ , 

/ 2 . 2* 
/x + y 

-yz/r /x 2 + y 2 


z/r 


0 

/ 2 , 2" 

Y x + y /r 


Performing the matrix multiplication gives the Cartesian components of 
acceleration as 




, 2 



a 

3 

/* Jj R eq 

X 

(5 sin 2 © -1) 

J2x 

2 

4 

r 

r 







a 

3 

x 

/4 R eq 

Z 

(5 sin 2 © -1) 

J2y 

2 

4 

r 

r 




2 



a J2z 

3 

2 

A 1 J 2 R eq 

4 

r 

z 

r 

(5 sin 2 © -3) 


where sin 8 = z/r. Note that the acceleration components given above 

are only the perturbing accelerations without the two-body term. 

The final calculation is the rotation of this acceleration vector 
from an equatorial to an ecliptic representation, 

a. J2 (ecliptic) = E £ J2 (equatorial), 
where E is the well-known rotation matrix involving the obliquity 


of the ecliptic 
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Now that all components of the perturbing acceleration, Sr. > 
defined, we can obtain £ r. and £ £. • To do this, £ r_ is numerically 
integrated with a generalized 4th Order Runge-Kutta algorithm for first 
order differential equations (Appendix 2). We can express $ £ as a 
set of first order equations 



£x can be numerically integrated to give £x(t) 
given the following initial conditions: 


Sri 

Si J » hen 


at t = t , 
o 



The propagation of ^r. (and ) continues until SjL i-S grester 

than or equal to some prescribed value Sr , then "rectification" 

max 

occurs . 

Hence, when £r ^£r at some time, t, 

y • mu v 


we reinitialize 




and compute a new reference conic orbit 0 Rectification ensures that 
the conic is always "close" to the true state* The propagation 
continues until £ r is again greater than or equal to £ r max > 
and so on to the end of the trajectory. 
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The frequency of rectification can be controlled by user input variable 
DRMAX. 

The use of constant thrust controls over discrete time intervals 
makes the trajectory discontinuous in acceleration at the control 
switching boundaries. The integration is thus done piecewise. That 
is, the state at the boundary between segments is used as initial 
conditions (rectification) for the next segment. 

4.2 Trajectory Termination 

There are four options for terminating trajectories in TRAJ: 

(1) final time, (2) closest approach to a target body, (3) sphere of influ 
ence the primary body, and (4) a radius relative to the primary 
body. Termination at final time is straight forward. For the other 
conditions, termination criteria are tested after each integration 
step. Once a cutoff condition is sensed, a step-size is computed so 
that TRAJ can propagate to an interpolated time. TRAJ takes the three 
previous planet relative position magnitudes plus the present relative 
position magnitude, and the corresponding trajectory times, and fits 
a third order polynominal through the four data points, using Newton's 
3rd order divided difference interpolation polynomial (Appendix 3). 

The independent variables for sphere of influence and stopping radius 
termination are the position magnitudes and the corresponding trajec- 
tory times are the dependent variables. Since the radius of the sphere 
of influence or the stopping radius is known before hand, the infor- 
mation that is needed is the time at these position magnitudes. For 
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closest approach, the same information is stored as before, but now 
the independent variables are the trajectory times and the position 
magnitudes are the dependent variables. Instead of knowing a time 
for which the position magnitude is a minimum, a value is computed 
corresponding to the minimum of a 3rd order polynomial. 



Figure (4-3) □ Radius Stopping Conditions 


4.3 Trajectory Accuracy 

As with all problems that require numerical integration, some 
criteria must be used in determining the nominal integration stepsize. 
The stepsize algorithm used in TRAJ is empirical, and it meets the 



Page 24 has been deleted. 
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requirements of reasonable numerical results and computer run time. 

The algorithm is 

h = € « j f| % 

where h is the integration stepsize, f is the gravity gradient (See 
Appendix 4) and € is a user input scale factor. There are con- 
straints on h such that 

* 

h ^ 5 days 

for heliocentric trajectories and 

h £ 1 day 

for planetary trajectories. 

The effects of h on the state transition matrix are small com- 
pared to the spacecraft state. Mass and mass variation are also not 
strongly affected by h because they are affected primarily by the 
spacecraft heliocentric position. Therefore, a good choice of h 
ensures a satisfactory trajectory. 

4.4 Trajectory Repeatability 

A major goal in building TRAJ was the ability to reproduce the 
same trajectory in all modes, given the same initial and spacecraft 
related data. To accomplish this, all propagation that resulted in 
altering the nominal integration stepsize was decoupled from the 
event stepsize logic. That is, event times and print times do not 
alter the norminal stepsize; instead, the information at the previous 
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nominal stepsize, t, is saved and a separate stepsize is computed 
for the event or print time, t + (Figure 4-5). If there are 
many events or prints between t and t + h a stepsize h„ is computed 
for each* Afterward, TRAJ returns to the nominal trajectory and 
continues the propagation. This suggests that, if the trajectory 
starts at any nominal integration step the trajectory will be 
duplicated, especially with respect to terminal conditions, for any 
run with different print and event times. For trajectories that do 
not start at a nominal integration step, there will be slight devi- 
ations in the terminal conditions from the nominal, but they will be 
close* 



^ - nominal integration time 

# - event or ‘print time 

event or print integration step 

h - nominal integration step 

Figure 4-5, Trajectory Preservation 
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4.5 State Transition Matrix Generation 

In a linear error analysis, the state transition matrix is used 
to map perturbations, 6x , about a reference trajectory at one epoch 
to perturbations at another epoch. First, the equations of motion 

* V 

x ” f (x, t) 

are expanded in a Taylor series around £ x , 

f (x + ix, t) = f (x, t) + S x + 0 (£x 2 ) 

9 

By neglecting powers of 0(£x ) and higher, the linearized differential 
equations 

$x = F £ x (4-1) 

are obtained, where 

$x = ^ (x + £ x, t) - f (x, t) 



The term f is a 3 x 3 time varying expression evaluated along the 
reference trajectory. The solution to the linear differential equa- 
tion 4-1, is of the form 
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£ is the gravitational acceleration contributon from secondary 
bodies, and ^ is the thrust caused acceleration term discussed in 
Section 4.1. % is a term of miscellaneous accelerations due to 

planetary oblateness, radiation pressure, etc. Since many of these 
accelerations have low order effects for Earth orbital missions, 

MAPSEP includes them according to user option via input control 
flags (See User's Manual.) It should be noted that all terms 
depend on the S/C position vector and this dependence must be 
taken into account when the variational partials are computed. 

To derive the variational equations, the vector function, f_, is 
expanded in a Taylor series about some reference solution to equation 
4-1. Hence, the right hand side of 4-1 becomes 

f (x+ Sx, t) = t (x, t) + ^ £ S x + & ( & x 

£ £ 

where <J X is relative to the reference trajectory state at time t. 

By neglecting second and higher order terms, this expression reduces to 

t 

Sx = F Sx (4-2) 

where is defined by 

&x = f (x + $x, t) • f (x, t) 


and F is a matrix of first order partials of the vector function f_ with 
respect to state components. Explicitly writing the F matrix, it is 
seen that it has only two non-zero partitions. 


F 


— 

1 

— 

0 

J 

I 

, _33 

1 

33 

f 33 

1 

°33 


1 



(6x6) 
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where is a 3x3 identity matrix and is a 3x3 matrix of time 
varying expressions evaluated along the reference trajectory. The 
components of f-j 3 > in terms of state variables, are obtained by 
analytically differentiating and 1 with respect to r. 

The solution of 4-2 is known from the theory of linear dif- 
ferential equations to be of the form 

Sx “ | Sxj, < 4_3 ) 

where $ is identified as the state transition matrix and is 
representable as 


4 


X, y, z,v x> v y> v z ) 

^ V V V V xo’ v yo, v zo^ (6x6) 


As noted before. & x is a perturbation, or deviation, from the 

-no 

reference trajectory at the initial epoch and, as such, it is arbitrary 
but constant. Differenting 4-3 with respect to time and substituting 
the resulting expression into 4-2 yields 
¥ 

f = F j (4-4) 

This equation is the variational differential equation for the state 
transition matrix and is numerically integrated to obtain $ over an 
arbitrary interval, t Q to t^, with initial conditions at t Q being given 
by 



From the more general point of view, equation 4-1 can be 



considered to be dependent on other parameters as well as the usual 
state variables. For example, the trajectory generated as a solution 
to 4-1 is dependent on thrust controls, gravitational constants of 
the Earth and sun, the value of J2 in the potential function, etc. 
When some of these parameters are to be investigated in a Linear 
analysis, the state vector of dynamic variables (x, y, z, x, y and 
is augmented to include the parameters of interest. This increases 
the dimension of the state vector to as great as twelve when there 
are three thrust controls, two gravitational constants, and J2 in 
addition to the vehicle state. 

For the problem where the state vector is augmented with para- 
meters, the equations of motion are more suitably written as 

-i* = i A (x, u, J2 , Jx P> JK g ; t) (4-5) 

where x ^ is the augmented state vector, i.e.. 


N* 
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In the above expression, x corresponds, as before, to the S/C state 
vector; u_ , to the thrust control vector; J2, the harmonic 

coefficient, and the f* ' s refer to gravitational constants. 

Following the previous analysis, it is possible to expand 4-5 to obtain 
variational equations for the augmented state transition matrix. This 
differential equation is of the form 

t* - f a (4 - 6) 

where ^ is partitioned as 
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Explicit definitions for the terms in are given in Appendix 4. Terms 
appearing in ^ are explicitly defined as follows: 


$ 


(partials of state component w.r.t. state components) 


4*0 


66 


© 

u 


(partials of state components w.r.t. thrust controls) 


4 H 


63 


0j2 = (partials of state components w.r.t. J2) 


3 £ 

TJ2 


61 


M = 
P 


(partials of state components w.r.t. the Earth* s . gravitational 
constant) = 



I 


61 


M - (partials of state components w.r.t. solar gravitational constant) 
s 

61 

Before concluding this section, it should be noted that MAPSEP 
has program logic which allows arbitrary augmentation of dynamic para- 
meters. That is, the program organizes and integrates matrices in 
equation 4-6 dimensioned to accommodate only those parameters requested 
during input* In this way, there is no time wasted in unnecessary 


calculations 
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4*6 Covariance Integration 

In any linear error analysis, a major problem is the prop- 
agation of state error covariances (P) from one event to the 
next event* Two methods are generally used: propagation with 

state transition matrices and numerical integration of the 
covariance matrix differential equations; both of which are 
options in TRAJ* In the latter case, the covariance is inte- 
grated to an event, where operations are performed on P, and 
the updated P is integrated to the next event* 

Given the nonlinear equations of motion 

X = X (x, u, <£> (4-7) 


where x is the spacecraft position and velocity, u are constant 
spacecraft controls and CU are time-varying thrust parameters 
(nominally zero), these equations can be linearized about a 
reference trajectory such that 


si = _Ai Sx + AA su + Ai 

^ X (J u 


Sil (4-8) 



31 


where £x, ^x, $u and <$6J are errors in the respective dynamic para- 
meters. Both and &c*?are described in terms of the 3x1 parameter 
set: pitch and yaw. The 6x1 ^ actually models two distinct 

processes for each parameter set (See also Page 62-A) . Whereas Equation 
4-7 describes motion of the deterministic reference trajectory, Equation 
4-8 describes the linearized propagation of trajectory deviations result- 
ing from dynamic and priori uncertainties. The covariance integrated 

by TRAJ not only maps dynamic errors but also measurement related 
errors, specifically, uncertainties in three station locations. The 
augmented state covariance is defined as 


p = e [ Si A £ s A ] 


where 



Sr 

$v 




^—'2 
S -3 . 


so that P = FP + PF + Q 

where F is similar to that used in the definition of the state transi- 
tion matrix, and is evaluated along the reference trajectory, and Q is a 
process noise matrix* The augmented F matrix is defined as 
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F * 


0 

f 

0 

0 

0 

0 

0 


I 

0 

0 

0 

0 

0 

0 


0 

g 

0 

0 

0 

0 

0 


0 

n 

0 

h 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 


where I is a 3 x 3 identity matrix. 


f - ^ 


dv 

dr 


S 


dv 


n = [gig] 


and h is the matrix of process noise correlation times 



Analytical equations for terms in the F matrix appear in Appendix A. 
The process noise, Q, is modeled as a stationary first order Gauss- 
Markov process and is defined as 
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0 0 0 

0 0 0 

0 0 0 

Q - 0 0 0 

0 0 0 

0 0 0 

0 0 0 


0 0 0 0 

0 0 0 o 

0 0 0 0 

-2 • h • E [ S W ] 0 0 0 

0 0 0 0 

0 0 0 0 

0 0 0 0 


For a discussion of Q, see Chapter 6 (GODSEP) • 

« 

Propagating P by integrating P is more mathematically 
accurate than the use of effective process noise as in the J method, 
and it lends itself to greater modeling flexibility for Q, The draw- 
back to this method is the increased run time. 
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5.0 TRAJECTORY GENERATION - TOPSEP 

The targeting and optimization mode, TOPSEP, generates a reference 
trajectory which is supplied as basic input to the error analysis and 
simulation modes. The primary purpose of TOPSEP is to incorporate in 
this trajectory all of the desired flight characteristics for a particular 
near-Earth mission while optimizing the final spacecraft mass. Initial 
conditions, a thrusting time history, and other control parameters are 
found which accomplish this optimization and yet lead to the required 
target conditions. The target constraints may be the final spacecraft 
state (cartesian or B-plane coordinates), final orbital elements, 
periapsis conditions, or other mission specifications which are listed 
in Table 5-1 and in the input section of the Users Manual (Volume II, P.18-19) 


Control Parameters 

Target Parameters 

Performance 

Parameter 

Initial State and Mass 

Geocentric State (Ecliptic 
or Equatorial) 

Final Mass 
(Payload) 

Thrust Magnitude 

Periapsis or Apoapsis 
Conditions (Radius , 
Velocity, Time) 


Thrust Times 

Orbital Elements 


Base Power Level 

Equatorial Latitude and 
Longitude 


Exhaust Velocity 

Final Mass 

i 



Table 5-1. Control, Target, and Performance Parameters 

The manipulation of trajectories to satisfy mission requirements is 
managed in three submodes of TOPSEP which represent successive stages of 
trajectory development. These submodes are: 
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1. Nominal trajectory propagation 

2. Grid generation 

3. Targeting and Optimization 

a. Trajectory targeting only 

b. A combination of trajectory targeting and optimization 

c. Trajectory optimization only 

Generally, these submodes are employed in order as listed above. However, 
any submode may be skipped or used individually if the proper control pro- 
file is available. Due to the simplicity of the first two submodes a 
brief discussion of their operational procedures is all that is necessary to 
understand their analytical basis in TOPSEP. The targeting and optimiza- 
tion submode will be reviewed in greater depth. 

5.1 NOMINAL TRAJECTORY PROPAGATION 

The simplest TOPSEP application is propagation of a single trajectory 
for spacecraft ephemeris information. After all the trajectory parameters 
are initialized, the trajectory is propagated from the initial state to 
the termination condition. TOPSEP performs no additional analysis of 
the trajectory when operating in this submode. This submode is also 
used for manual manipulation of the control profile. 

5.2 GRID GENERATION 

The grid generation subraode is available to produce a number of tra- 
jectories which do not necessarily satisfy mission requirements but pro- 
vide a range of trajectory solutions. Thus, the main purpose of the grid 
submode is to locate desirable control regions for further examination. 
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In turn, each control Is incremented a fixed amount while the remaining 
controls maintain their nominal values. A single low thrust trajectory 
is generated for each control change and the associated target error index 
is calculated. Subsequently, contours of constant target error may be 
plotted in the control space so that some control regions can be eliminated 
from further consideration. Upon completion of the grid, the trajectory 
generation mode is terminated and the program user must choose the best 
control profile to initialize targeting and optimization or to employ 
another grid approach. 

5.3 TARGETING AND OPTIMIZATION 

The trajectory targeting and optimization submode features a discrete 
parameter iteration a Igor ithm which accommodates the non-linear aspects of 
the low thrust problem. The algorithm is a modification of Rosen’s pro- 
jected gradient method (PGM) for non-linear programming (Refs. 5 and 6). 

The parameters which have been chosen to shape the trajectory (Table 5-1) 
constitute the control profile and are subject to modification by the PGM 
algorithm. Based upon the sensitivities of the final S/C mass and state 
to control variations, corrections to the profile are computed which 

maximize performance while minimizing target errors. The performance is 

o 

measured simply by the value of the final spacecraft mass while the target 
errors are measured according to the constraint violations. The method 
chosen to represent the target errors in terms of a scalar measure is the 
quadratic error index which is the weighted sum of the squares of the 
target errors. 

When the targeting and optimization submode is entered, a nominal 
trajectory is propagated directly from the input parameters. A series of 
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tests is performed to determine which submode (targeting, optimization or 
both) is to be executed. If the target error index is large, the submode 
will be exclusively targeting. However, a target error index smaller than 
some value (TUP in namelist $T0PSEP) will result in simultaneous targeting 
and optimization. Whenever the index is below a specified lower bound 
(TL0W in namelist $T0PSEP), the optimization algorithm will be executed. 

Prior to the application of the projected gradient algorithm, the 
targeting sensitivity matrix S and the performance gradient are computed. 
Elements of the S matrix represent the sensitivities of individual target 
parameters to changes in controls and are used for both targeting and 
optimization. Similarly, the elements of the £ vector represent the 
sensitivity of the performance index to changes in controls although 
these elements are used only for optimization. For purposes of targeting 
only, S is computed from the integrated state transition matrix (STM) 
or by finite differencing techniques and ^ is ignored. Appendix 7 
discusses the formulation of S from the integrated STM. Whenever opti- 
mization is to occur, both S and £ are constructed by finite differenc- 
ing techniques. Following' the determination of S and a weighting 
matrix which amplifies or diminishes the effects of the chosen controls 
is calculated. Applying the projected gradient algorithm a control cor- 
rection is established. The magnitude of the control change is determined 
by computing trial trajectories. The new control profile is simply the 
old control profile plus a scalar multiple of the control correction 
such that the targeting error index is minimized and/or the performance 
index is maximized. If the optimization is complete (the values of the 
performance index have converged to a maximum), TOPSEP is terminated. 
Otherwise, the submode decision is made again and the cycle is repeated. 


* 
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Of primary importance in the targeting and optimization submode is 
the selection of the control correction. In the following sections this 
selection process will be discussed. 

THE PROJECTED GRADIENT METHOD 

The projected gradient method has been devised to maximize a per- 
formance index while simultaneously minimizing an error index. Since 
maximizing performance is equivalent to minimizing cost in an optimization 
sense, PGM 1 s purpose relative to the trajectory problem may be restated as 
minimizing fuel expended as well as minimizing target error. The concept 
of net cost, which is simply a more realistic assessment of fuel expended, 
will be discussed later in this section. 

The projected gradient algorithm employs cost-function and constraint 
gradient information to replace the multi-dimensional targeting and opti- 
mization problem by an equivalent sequence of one-dimensional searches 
(Ref. 7). In this manner, it solves a difficult multi-dimensional problem 
by solving a sequence of simpler problems. In general, at the initiation 
of the iteration sequence, PGM primarily satisfies the constraint require- 
ments. As the iteration process proceeds, the emphasis changes fron con- 
straint satisfaction to cost-function reduction. 

Since numerous analytical developments of this technique are available 
(Refs. 5 and 6), this presentation will primarily emphasize the geome- 
trical aspects of the algorithm. Clearly, the geometric interpretation of 
the algorithm is the motivation for the logic contained in TOPSEP, and a 
basic understanding of these concepts is usually sufficient to enable the 
user to efficiently manipulate ^TOPSEP input and to handle diverse mission 


problems . 
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PROBLEM FORMULATION 

The projected gradient method solves the following non-linear pro- 
gramming problem: 

Determine the values of the control variables, ij, that minimize the cost 
function (optimization variable) 


F(u) 

subject to the equality constraints 

e (x(u)) = T (x (u) ) - T d = 0, 

u £ U, an M-dimensional control space 
T £ T, an N-dimensional target space 
x £ X, a six-dimensional state space 
M > N 

where the elements of e, T, T,, and x are referred to as the target error, 

~"U 

the target values, the desired target values, and the final state respec- 
tively; and F is a scalar valued function measuring system cost. 

In an attempt to solve the constrained optimization problem, iterative 
methods are employed. The following scheme briefly describes the process 
occurring in TOPSEP. 

• Guess (In general u^ will not satisfy the constraints nor 
minimize F) . 



AO 


• Determine a correction £ 11 such that 
•• F(u q + 4 u) <. F (u^) and 

l! +A!i)|| < || e (u^)] 


• • 


• Iterate until 

•a F is minimized and 

•• I! —I ! 44 £., a pre-determined tolerance 
NOMENCLATURE AND CONCEPTS 

To fascilitate the discussion of the projected gradient algorithm, 

the following nomenclature and basic concepts will be introduced, x denotes 

a column vector whose elements are x^, where i = 1,2, ••*, n and n is the 

T 

dimension of the space containing x- Y denotes the transpose of the real 
matrix The feasible region defined in the M-dimensional control space 

within which PGM operates is the restricted space 


Vn - - u ‘ 


MAX 


, i * 1, 2, •••, M. 


The equality condition implies that the control is on a bound. The cost 

gradient g is an M-vector of partial derivatives and is defined as 

T 




The sensitivity matrix is that matrix whose rows are the gradients to the 
equality constraints, and is denoted by 

S(u) = ±t^L. 

d a ■ 

where e is an N-dimensional vector. The target error function is defined 
to be 

E<u) = e T [wj e 
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where W is a target weighting matrix which will be defined later, 
e 

Corresponding to each control vector u in the control space U, 
there is an error vector e_. Let A q be the set of all u such that 

e(u) “ 0. 

A then represents all the control vectors satisfying zero-target error, 
o 

It can be shown (Reference 6) that A q defines an M-N dimensional non-linear 
hypersurface or manifold in U. Unfortunately, A q cannot be defined 
explicitly; hence, one cannot easily find a u which is an element of A q . 

However A can be estimated implicity via the sensitivity matrix. 

9 o 

Let A be the set of all u such that 
c 

e(u) = c, 

where c is a vector of constants. Thus, A c represents a non-linear mani- 
fold containing those control vectors which provide constant target error. 
It can also be shown (Reference 6) that any u in the control space is 
contained in one and only one A^. At a given u, the corresponding 
non-linear manifold A £ may be approximated by a linear manifold B(u) which 
is defined explicitly by the sensitivity matrix S(u). The linear manifold 
B(u) may be considered a tangent hyperplane to A c at u. The orientation 
of B (u) in the control space allows one to define a search direction to A q 
which is orthogonal to B(u). This search is in the direction of maximum 
decreasing target error • 

Let B(u) denote the orthogonal complement to B(u). One can demon- 
strate (Reference 6) that B(u) is the unique linear space that can be 
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translated to obtain the linear manifold B (u) . Furthermore, there 
exist unique orthogonal projection operators P(u) and P(u^ that resolve 
any vector in the control space into its corresponding components in 
B(u) and B(ti), respectively; that is 

u = P(u) u + P(u) u. 

In particular, 

P(u) = S T (SS T ) S (5-1) 

and 

P(u) = I - 'P 

A* 

where I is an identity matrix. The pro jec tion operators p and P thus pro- 
vide a simple method for reconstructing a general vector Au emanating from 
u into its two components in B and B. In a discussion to follow later in 
this section, it will be explained how Au may be defined such that PAu 
represents the control correction to minimize the target error and PAu 
represents the control correction to minimize cost. 

The final key concept employed by PGM is the idea of problem scaling. 
The purpose of problem scaling is to increase the efficiency of the tar- 
geting and optimization algorithms by transforming the original problem 
into an equivalent problem that is numerically easier to solve. 

To numerically scale a problem, two general types of scaling are 
required: 1) control variable scaling, and 2) target variable scaling. 

Control variable scaling is accomplished by defining a positive diagonal 
scaling matrix, (UWATE in namelist $T0PSEP), such that the weighted 

control variables are given by 

u f = W u . 

— U “ 
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Similarly, target variable weighting is accomplished by defining a 
positive diagonal scaling matrix, (IART0L in namelist $T0PSEP) , such 
that the weighted target variables are 

T' (u) = [wj T (W u -1 u') . 

The target error index is then 

E(u) = e' T e' . 


TOPSEP contains several options for computing the control variable 
weighting matrix. The option most often used is the normalization scaling 
matrix (See Appendix 6 for other options). 



The target variable weighting matrix is always computed as the reciprocal 
of the constraint tolerances and is given by 


w 


e i ii 


£ i 


th 


where is the tolerance for the i target error. 

For simplicity, the following discussion of the algorithm assumes an 
appropriately scaled problem. However, the scaled equations can be obtained 
by making the following simple substitutions. 

u replaced by u 1 
e_ replaced by e l 


S replaced by 
replaced by 


Ou ]' 1 * 
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DIRECTION OF SEARCH 

The concept of the direction of search in control space needs 
slightly more elaboration. The direction of search is nothing more than 
a particular line in the control space along which the target error is 
reduced or along which the cost function is decreased. In a more precise 
sense, the direction of search at u is a half-ray emanating from u. Thus, 
for any positive scalar, if , the equation 

A \/ ft 

u = u + OAu 

sets the limits of this half-ray and represents a "step" in the direction 
Au from u . This is illustrated in Figure 5-1. 



This concept of direction-of-search is particularly important since it 
enables the M-dimensional non-linear programming problem to be replaced by 
a sequence of one-dimensional minimizations. What remains to be explained 
is: 1) how to select the direction of search, and 2) how to determine 

the step size in that direction. 
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The projected gradient method uses two basic search directions. 

For this discussion they will be termed constraint and optimization 
directions. PGM proceeds by taking successive steps in one or the other 
of these two directions. The computation of each of these search direc- 
tions is described below at a particular point « in the M-dimensional con- 
trol space where N constraints (target conditions) are enforced. 

constraint direction 

Th. constraint direction depends critically on the nnmber of targets. 

Two cases are distinguished below: 

1. If N<M, then that unique control correction A u is sought which 

solves the linearized constraint equation 

[s]au + e (u ) = 0 < 5_2) 

and minimizes the norm of -*u . The -solutions to the preceding 
vector equation define the M-N dimensional linear manifold B q (u) 
which is an estimate of the non-linear manifold A q (zero-target 
error). The desired minimum norm correction A u is then the vector 
of minimum length in the control space from u to the linear mani- 
fold B (u), thus requiring 4u to be orthogonal to B q (u). 
o 

Application of the linear operators P and P allow one to represent 

A 

Au as the sum of two orthogonal vectors relative to B q (u) or 

iiu = P Au + P ^ U. > 

however > 

P Au - 0 

since there are no components of in B q * From equation (5 X) 

\x may be expressed as 
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* T T a 

A u = S (SS ) S4u 

which may be reduced to 

du = -S T (SS T ) e(u) 

using the constraint equation (5-2). This correction is illus- 
trated in Figure 5-2. 



Figure 5-2. Illustration of Minimum Norm Constraint Direction 
for N = 2, M * 3. 

The direction of search then is simply taken to be this minimum 
norm correction to the linearized constraints, 

2 # If N = M there is unique solution to the linearized constraint 
equations without the additonal requirement that the norm of 
the control correction be minimized. The solution for Au reduces 
to the familiar Newton-Raphson formula for solving M equations 


with M unknowns; namely 
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Au. = " S 1 e(u) 

The Newton-Raphson correction is illustrated geometrically in 
Figure 5-3. 

Second linearized 



Figure 5-3. Illustration of Newton-Raphson constraint direction for 
N = M = 3 

OPTIMIZATION DIRECTION 

When the number of targets is less than the number of controls, it 
is then possible to minimize the cost function F (u ) assuming, of course 
that u is some non-optimal control profile. Obviously, the steepest 
descent direction, - g <u) , would be the best local search direction for 
reducing the cost function. Such a direction, however, could produce 
unacceptable constraint violations. To avoid this difficulty PGM ortho- 
gonally projects the unconstrained negative gradient,- £, onto the local 
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linearized constraint manifold B (u) . By searching in the direction of 

c 

this negative projected gradient the algorithm can guarantee in a linear 

A 

sense that there is no further constraint violation than that of e(u) . 

To calculate this direction, it is only necessary to apply to - g the 
projection operator P(u) which will map the vector into its component on 
the linear manifold B^ (u) • Thus, 

A A 

Au = - P g (u) 

-l - 

= -[l - S T (SS T ) SJ g (u) 

COMBINED TARGETING AND OPTIMIZATION DIRECTION 

When it is desirable to minimize the cost function as well as reducing 
the target error the constraint direction and optimization direction may 
be combined such that the resulting control correction is of the form 

Au = + A u^ 

where Au^ is the optimization correction and is the constraint correc- 

tion. Note that A andAUj are orthogonal components of Au . Depending 
upon the magnitude of the target error, one may want to emphasize either 
optimization or targeting. Since this decision is rather subjective and 
linked directly to the degree of problem non-linearity, an option is pro- 
vided to weight Au (See Page 50) in one of its component directions. 

Figure 5-4 and Figure 5-5 illustrate the geometric interpretation of the 
resulting control correction. 



Nonlinear Constraint 
Manifold (e=0) 


Figure 5-4. 


Constrained 

Optimum 


Geometric Interpretation of a Combined Targeting 
and Optimization Control Correction, N-l, M-3 



B (u, 
c — 


B (fi) 
c 


Figure 5-5. 


Illustration of Combined Targeting and Optimization 

Control Correction As Seen In B (u), 

c 

y A 

(Enlargement of B c (u) from Figure 5-4). 
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The total control correction is constructed as follows, where d is an input 

scalar ' A » = IMI * a * 

il*% II 

(DP2 in namelist $T0PSEP) * Thus , one has the flexibility of determining 
the magnitude of the optimization correction relative to the magnitude of 
the constraint correction* The optimization correction can then be written 



The norm of the control correction ^u^which is obtained by summing^ii^ and 
A not as important as the direction* The resulting half-ray provides 

the basic search direction with which to calculate the trial step. 

TRIAL STEP-SIZE CALCULATION 

At any particular point u in the control space, the PGM algorithm 
proceeds by reducing the multi-dimensional problem to a one -dimensional 
search in the direction prescribed by the constraint or optimization con- 
trol change vector. Once the initial point u and the direction of search 

A 

Au are specified, the problem reduces to the numerical minimization of a 
function of a single variable;, namely , the step scale factor^ . PGM per- 
forms this numerical minimization by polynomial interpolation based on 
function values along the search ray and the function's value and slope at 
the starting point. 

CONSTRAINT DIRECTION 

The function to be minimized along the constraint direction Au^ is 
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E(tf), the sum of the squares of the target errors. 


) = e T (u+Au) tO e 


Evaluation of the function at ? = 0 results in 


E (0) = ej (u) w e (u) 


Differentiation via the chain rule yields 


3e(* 


T A A 

= 2e (u) S A u 2 


o 


If constraints are reasonably linear, a good initial estimate for the 
minimizing is one. 


OPTIMIZATION DIRECTION 


The function to be minimized along the optimization direction is 
the estimated net cost function 0(3?), where 


C(« ) *= F(u + SAuj) - F(u) + 


it- 


S T (SS T ) e (G + $ A U^ 


Change in cost produced 
by a step of length 
i ||AUi|i along 


Linearized approximation to change 
in cost required to perform a minimum 
norm correction in order to retarget 


Clearly, 


c(0) = - g S T (SS T ) e (u) 


By expanding C(0) in a taylor series in *8 about ^ = 0, and by making use 
of the fact that P = 0 , it can be shown that 


= g A u 


% - 0 


These properties are illustrated in Figure 5-6. 
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Both the constraint and optimization directions are based on a 
sensitivity matrix assuming a linear space. Due to nonlinearities in the 
problem, it is often necessary to restrict the trial step such that 
unexpected increases in cost or target error are reduced. In addition, 
the control correction Au does not reflect the "nearness" of control bounds 
as long as u is not on any bound. Thus, the trial step must also be 
restricted so that the new control vector remains within the bounded con- 
trol space. For these reasons a maximum limit is placed on "6 . After 



Figure 5-6. Properties of net-cost function along the direction of 
search. 
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i s computed, a calculation is made to find the along the 
search direction which will size a trial step that intersects the 
boundary of the feasible region. This value is compared to the 
maximum step allowed by the user to counter the nonlinearity prob- 
lem. The smaller value is specified as the maximum 'ZT allowed 
during the search. A method has been devised to alleviate the 
problem of a control which is very near a boundary. A tolerance 
region is defined in the neighborhood of the boundary surface such 
that if the control is within this region and A u intersects the 
boundary the search to minimize the net cost or target error can 
continue along the boundary without calculating a new sensitivity 
matrix. Once a control element reaches one of its bounds it becomes 
inactive. Unless a subsequent correction for this control element 
is back into the feasible region it remains inactive. 

ONE DIMENSIONAL MINIMIZATION 

Nonvariant minimization in PGM is performed exclusively by 
polynominal interpolation. The actual function to be minimized is 
fitted with one or more quadratic or cubic polynominals until a 
sufficiently accurate curve fit is obtained. -The minimum of this 
curve and the corresponding scale factor can easily be found analyt- 
ically. 

The one-dimensional search proceeds by taking trial steps in 
the Au direction to obtain information about the function to be 
minimized. If VAu is a constraint correction^ the quadratic 
error function is evaluated; if TfA u is an optimization correction , 



54 


the net-cost function is evaluated; and if TAu is a combined cor- 
rection, a function which is a weighted combination of the error 
function and net-cost function is evaluated. 

The minimization routine makes ingenious use of all the infor- 
mation it accumulates. The following curve-fitting techniques are 
applied in order. 

1. Quadratic polynomial fit: two points-one slope; 

» 

2. Cubic polynomial fit: three points-one slope; 

3. Quadratic polynomial fit: three points; 

4. Cubic polynomial fit: four points. 

Each time a trial step is taken, the function which is evaluated is 
used as a trial point to analytically determine the next trial step. 
The analytical formulation of the above curve fits may be found in 
the subroutine description of MINMUM in the Program Manual . 
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6. LINEAR ERROR ANALYSIS - GODSEP 

GODSEP analyzes spacecraft and trajectory related dispersions as 
a function of expected uncertainties in dynamic and navigation parameters. 
The ensemble of expected errors is studied without actually simulating 
individual trajectories by applying linear techniques. That is, small 
deviations about a reference trajectory are linearly related to other 
deviations by a transformation matrix. For example, the state transition 
matrix relates position and velocity deviations about the reference tra- 
jectory from one time point to another. The ensemble of errors, or covar- 
iance, is assumed to have a zero-mean Gaussian distribution, except for 
special processes. 

Probabalistic a-priori errors in the environment, spacecraft and 
tracking systems are propagated in time along the reference trajectory 
through sequential events such as orbit determination (OD) and guidance 
corrections. Two types of ensemble error or covariances are distinguished 
knowledge , which reflects the ability of the OD algorithm to estimate the 
spacecraft state and other parameters; and control , which represents the 
dispersions of the actual spacecraft trajectory about the reference. 
Covariance propagation is done by either integration of covariance varia- 
tional equations, or by the state transition matrix method. 

The error analysis proceeds sequentially from start time through each 
specified trajectory event to final time. Event types available are 
measurement, propagation, eigenvector, prediction, thrust switching, and 
guidance. A measurement event processes tracking data at a time point by 
applying the user specified OD algorithm. Available to the user are both 
Kalman -Schmidt (K-S) and sequential weighted least squares (WLS) filters. 
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GODSEP modularity also allows the user to insert his own filter algorithm 
quite easily. The filters are distinguished by their methods of gain 
matrix calculation and subsequent update of the knowledge covariance. 

A propagation event merely updates the knowledge (and control) covari- 
ance at the event time. Its primary value is in maintaining accurate covari- 
ance values during long propagations by forcing computation of the effective 
process noise over predetermined, user-specified intervals. 

An eigenvector event is used for information display and behaves similar 
to a propagation event. Covariance matrix sub-blocks are converted to 
standard deviations and correlation coefficients. It also computes eigen- 
values, their square roots, and eigenvectors for the position and velocity 
3x3 sub-blocks of the state covariance matrix. Thrust switching events are 
simply eigenvector events at the time where a change in the number of 
thrusters or thrust policy has occurred. 

A guidance event is an update of the control covariance to reflect 
implementation of a trajectory correction. A correction is not performed 
deterministically, but only in a probablistic sense. The guidance event 
computes expected correction covariances (Av or thrust control), target 
error covariances before and after the guidance event, and the updated state 
control covariance. 

The following sections will describe in more detail the analytical 

foundations of GODSEP. 

6.1 AUGMENTED STATE 

The augmented state discussed previously in TRAJ (Section 4.5) 
includes dynamic parameters besides the basic spacecraft position and 

In GODSEP, the augmentation not only adds measurement 


velocity vectors. 



57 


related parameters to this list, but also distinguishes between solve- 
for and consider. Solve-for parameters are directly estimated by the OD 
process* Consider parameters are system uncertainties which are recog- 
nized and accounted for in the estimation algorithm but are not estimated, 
usually because the process cannot be adequately modeled or there is a high 
correlation between two (or more) parameters which might cause numerical 
difficulties if both were solved-for. 

The possible augmented parameters that can be either solved-for or 
considered are 

r • thrust bias (magnitude and pointing) 

dynamic - • J2 zonal harmonic 

, • gravitational constants of primary body and/or sun 

measure- f • tracking station locations 
ment ] 

l • sensor bias (range, range-rate, etc.) 

Time varying thrust noise (magnitude and pointing) can only be considered 
in the standard GODSEP analysis, but can be solved-for (or considered) in 
the covariance integration option (PDOT in Section 6.2), A third possible 
category, in addition to solve-for and consider, is the ignore parameter 
used in generalized covariance analysis (Section 6.5). 

The total ensemble of state uncertainties, or error covariance, 
including all augmented parameters, is formed by applying the expected 
value operation on state deviations from their reference values, 

P = E £ $x S x T J 

The covariance P contains uncertainties and their respective correlations 
for all parameters in the augmented state. There are two covariances, 
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corresponding to control and knowledge, which are computed in parallel 
through a GODSEP analysis. Starting with a-priori values, each P is 
modified between events by trajectory propagation effects, and at events 
by either OD (knowledge) or by guidance corrections (control). 

6.2 COVARIANCE PROPAGATION 

There are two methods. available in GODSEP for propagating covariances 
between events: transition matrices ($) and explicit covariance integration 

(PDOT) . Although these two techniques were discussed in TRAJ, Sections 4.5 
and 4.6, respectively, their importance in GODSEP requires additional explana- 
tion. 

The most common form of covariance propagation, both in GODSEP and in 
other linear error analyses, makes use of transition matrices. This is 
because the <f's are a characteristic of the trajectory, not of the covari- 
ance. A covariance P(t^) is propagated from time t^ to, t^ by 

P(t 2 ) - P(t x ) + Q 21 (6-D 

where * n = 4 (t^) and Q n = QOyt^) is an effective process noise 
covariance. 

Transition matrices can be stored, for example on magnetic tape, to 
be used for analyses of different error source levels and navigation strate- 
gies. Obviously, if *’s have been computed between the intervals t Q , t^ 

t •••t they can be used to propagate any P as long as the set of pro- 
2* N’ 

pagation times is a subset or equal to the original set. Transition matrices 
can always be chained to cover desired propagation intervals, for example, 

rf 3 l = «f(t 3 ,t x ) = *(t 3 ,t 2 ) *( t 2 ,t x ) 
then P(t 3 ) = 4 31 P(t t ) 0 31 + Q 31 


letting 
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The method of computing and storing 4 ' s (on the STM file) over a grid 
of time points is used in GODSEP to facilitate parametric error analyses. 

Since covariance propagation accounts for uncertainties in all 
dynamic parameters which have been augmented to the basic spacecraft state, 
the transition matrix must have the same augmentation. In actual operation, 
TRAJ provides a transition matrix containing only dynamic variations which 
GODSEP must augment with appropriate rows and columns of zeros (for measure- 
ment parameters) such that the total augmented 4 is consistent with P. 

An additional requirement for GODSEP is the modification of the thrust 
sensitivity matrix 0 computed by TRAJ as part of the augmented transition 
matrix. 


0 (t^) 


^S.(t 2 ) 

Th 


where u are constant thrust ’controls (proportionality, cone and clock) over 
the interval (t ,t 2 ). The 0 matrix is used in GODSEP to map thrust biases 
into spacecraft state uncertainties. However, GODSEP thrust biases refer 
to a single thruster. If more thrusters are operating, and if each operating 
thruster is assumed to be independent of all others in terms of bias, then 
the total effective bias is simply the single thruster bias divided by i/IT 
where N is the number of thrusters, or 



The effective process noise, Q , is a very important conditioning 
term on the covariance propagation. Because a rigorous mathematical com- 
putation of Q involves (1) modeling of a process or processes which are 
ill-defined and (2) evaluation of complex double integrals, GODSEP uses a 
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simple analytic approximation. The effective process noise assumes that 
time-varying thrust errors appear as stationary first-order Gauss Markov 
processes. The more rigorous modeling is performed in the PDOT option to 
be discussed shortly. The relationship between P and Q precludes the aug- 
mented state from containing time-varying thrust terms so that process noise 
takes on the appearance of consider parameters. 

The effective noise over a time interval t^ to t^ directly affects 
only the spacecraft position and velocity uncertainties at 


*21 


’/* M _ /“o ' |t 

|L v + 121 lo v V 


where J is the 6x6 state transition sub-matrix of the augmented transition 
matrix (^) , 

H x = g(t x ) P<* g T (t 1 ) 

= g(t„) g T (t„) 




T 1 < y /N 

0 

0 


0 0 

2 . 


t 2 c 2 /n o 
0 


T 3 a- 3 /N 


oL = 


ot 


0 

*0 


0 


S 2, if At>T i 
0, if At A T ± 


At 


_ 2 ~ 2 — 2 
^1 ’ ^2 1 ^3 


N = 
g = 


C 2 • h 


variances in thrust proportionality, cone, clock, 
respectively, 

number of operating thrusters 

(see Section 4,1 and Appendix 4) 


T l> T 2> T 3 


= process noise correlation times 1 
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Thrust proportionality error is scaled by the number of thrusters 
because it is assumed that time-varying noise is independent for each 
thruster, just as bias is. Thrust pointing noise is also scaled 
because it is assumed to be caused by the thrust vector control system 
which currently consists of gimbaling each thruster independently. 

This empirical model for Q is generally effective over propagation 
intervals on the order of 2 days or less. Propagation events can be 
employed in GODSEP to break up longer intervals and to ensure the accuracy 

of. Q. 

GODSEP has the capability to model two simultaneous processes 
in each of the three noise elements: thrust proportionality and 

two orthogonal angles (such as pitch and yaw). The first process is 
always assumed to be thruster dependent whereas the second process can 
be either thruster dependent or independent, A thruster dependent 
process is proportional to the number of operating thrusters, e.g., 
errors in thruster alignment and beam divergence. A thruster indepen- 
dent process occurs no matter how many thrusters are operating, e.g., 
errors related to the inertial reference system. 

For Earth orbital missions, a potentially important process is 
the thrust start-up sequence after exiting from solar occultation. 

A small but significant amount of time is required to restart the 
thrusters depending upon thermal conditions and thruster design. 

Although thruster delay time can be included during the reference 
trajectory generation (TOPSEP), there are still uncertainties in 
thrust initiation caused by errors in thermal modeling and 
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calibration, etc. Thrust restart error is modeled in GODSEP as an 
effective start-up time uncertainty. The effective process noise 
matrix, Q, at nominal thrust start time is 

2 T 

Q = E C* E 
H & on “■ 

where & is the sensitivity of vehicle start with respect to thrust 
on time, that is, g. O x/ and is the varlance of 

thrust on time, that is, O ' = E jT A t Qa t Qn J • 

The second method of covariance propagation, PDOT, is used primarily 
to examine thrust noise effects. Process modeling is mathematically 
rigorous and includes augmentation of thrust noise parameters to the basic 
state. Recalling from Section 4.1, the linearized equations of motion, 

£ x = F Sx 

and from Section 4.6, the corresponding covariance matrix differential 
equations 

P = FP + PF T + Q 

where x is the augmented state which, for PDOT, includes at most space- 
craft position and velocity, thrust biases, thrust noise and tracking sta- 
tion locations, F is the variation matrix, and Q is a white noise term 
which affects only the thrust noise directly. If thrust noise is omitted, 
then the integrated covariance would in theory be identical to a similarly 
augmented covariance propagated by transition matrices. 

In PDOT, the time-varying noise is modeled as a stationary Gauss- 
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Markov process, as in Q, 

oO = 

where cO is a 6x1 vector of independent noise parameters corresponding to 

thrust proportionality, cone and clock, each of which is described by two 

processes having their own distinct correlation times (T) . This permits 

the study of superimposed and multi-process effects. W is a white noise 

component which drives the time varying noise (and defines the only non- 

T 

zero term of Q which is E [W W ] ). 

Since CO is in the desired form of the linearized equations of motion 
it can easily be augmented to the state vector (and covariance). Thus, 
to can be solved-for in the PDOT mode although in reality this practice is 
questionable because of the to modeling assumptions - who knows how thrust 
noise really behaves? 

One of the more useful applications of PDOT is in refining the form 
of effective noise, Q, for a particular mission and in verifying the 
explicit assumption in Q of zero correlation between noise and state para 
meters . 

Whether state transition matrices or PDOT is used for covariance 
propagation, an auxilliary computation is the vehicle mass uncertainty. 
Since mass and thrust magnitude uncertainties are indistinguishable 
in their trajectory effects, that is, they are correlated one to one, 
GODSEP has chosen to model thrust (acceleration proportionality) magni- 
tude explicitly, and provide the approximate equivalent mass uncertainty 
as supplementary information. 


1 • 


T, 


CO + W 
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Two types of mass uncertainty are distinguished: knowledge 

(estimated) and control (actual). Estimated mass uncertainty is the 
instantaneous knowledge error in thrust magnitude, bias and noise. 

estimated (T 2 = ( (T . 2 + (T 2 ) 

m a D 3tl 

2 2 2 

where (T (7* , , (T are the variances in mass, thrust bias 

m ’ v ab an 

proportionality (from P ) and thrust noise proportionality (from P^ 
or Q), respectively. 

Actual mass uncertainty is the cumulative mass variation reflected 
by the control error covariance. The actual mass deviation from the 
reference at time t + A t based upon uncertainties from time t is 

actual <r m 2 (t + 4 o -[<r.(t>+! <r ab At] 2 + 2 <r sn 2 T 2 » 

* 1‘- --¥] [-rf 

where m, (T. and (T are averaged ever the interval At, c is the 
w ab an 

exhaust velocity and T the correlation time* (and for PDOT) 

are obtained from the augmented control covariance. Accuracy of the 

mass variance computation depends upon the event schedule because GODSEP 
_ 2 

evaluates \j from one event to the next. 
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6.3 MEASUREMENT TYPES 

In a linear error analysis, the reference trajectory deterministically 
characterizes the motion of the S/C, and no real state vector estimation 
is explicitly performed in GODSEP. Rather, an orbit determination anal- 
ysis estimates how well the state vector can be determined if the S/C 
were to move along the reference trajectory and were to be observed as 
directed by the analyst. In this sense, the term ’’orbit determination" 
refers to the calculation of a knowledge covariance based upon the proc- 
essing of modeled observational data. This section of the Analytic Manual 
describes the data types and mathematical models that have been implemented 
in GODSEP. The next section will treat the problem of filter formulation 
and the process of updating the knowledge covariance. 

When an observation is to be simulated in GODSEP, the knowledge 
covariance is propagated to the scheduled measurement point and is made 
available to be updated by the filter. Before this can happen, it is 
necessary to evaluate the observation matrix which relates the observables 
to the state vector. Given an arbitrary vector (or scalar) measurement 
X = i QD where X is t ^ xe tota l augmented state consisting of 


► < 

X 


* * 

spacecraft position and velocity 

s_ 


solve-for parameters 


- 

dynamic consider parameters 

v 


measurement consider parameters 

w . 

L — J 


ignore parameters 


then the linearized measurement, which assumes small deviations from the 
nominal, is 

Cy = H Sx + H £ s + H Su + H £ V + H £ w 
° X s — U V w 
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where H 25 
x 


^ z 


3 x w ^ w 

all of which are computed analytically in GODSEP. 

GODSEP has the capability of processing the following measurement 
types : 


V Z 


are the observation matrices 


earth- 

based 




L 


o 2-way range 
o 2-way doppler 
o 3-way range 
o 3-way doppler 

o simultaneous 2 -way and 3-way range 
o simultaneous 2-way and 3-way doppler 
o differenced 2-way and 3-way range 
o differenced 2-way and 3-way doppler 
o azimuth and elevation angles 


spacecraft 

based 


o horizon scanner angles 
o star-planet horizon angles 
o apparent planet diameter angle 


For Earth based data types, nine tracking stations are available. 
The default station locations are intended to provide global coverage 
with Fairbanks, Alaska having the highest geographic latitude (64°) 
and Canberra, Australia at the lowest latitude (-35 ). To reduce 
the amount of Earth based tracking, spacecraft based data types are 
often used in a semi- or fully autonomous system* 



Page 64-B has been deleted. 
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The following definitions of position and velocity vectors are neces 
sary to relate the mission geometry to the observable quantities. All 
vectors are assumed to be column vectors and are expressed relative to an 
inertial, ecliptic coordinate frame, unless otherwise noted. 


S/C heliocentric cartesian position and velocity 


Earth heliocentric cartesian position and velocity 


-l’-l 


Station 1 geocentric cartesian position and velocity 


- 2’ -2 


Station 2 geocentric cartesian position and velocity 


/O /0 
'-l’-l 


/O /© 

- 2 ’-2 


S/C position and velocity relative to Station 1 


S/C position and velocity relative to Station 2 


Pi 



Unit vectors defining direction of S/C from 
Stations 1 and 2 respectively 




Ap = 
-2 


S/C range and range-rate from Station 1 
S/C range and range-rate from Station 2 
3-way range and range-rate 

Differenced 2-way and 3-way range and range-rate 

Ceocentric spherical or cylindrical coordinates of 
Stations 1 & 2, s = (r, or (x > ) 

Zero vector, 3x1 

Identity matrix, 3x3 unless noted otherwise 



S/C 


Figure 6-1. 


Tracking Geometry for Range and Range-Rate 
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We first note one identity which is used numerous times in the follow- 
ing derivations. Given the vector a. = b. - c_ and its corresponding 
unit vector S = a/ | a. | > 


^a 


2 b j a | 


[ i - 


* a. T 1 

a a J 


( 6 - 2 . 1 ) 


Two-way range and range rate from Station 1 are modeled 


where 


rn rs 

<°i ■ £ Pi 


p - e T P 

'l L 1 1 


-1 


= x - 




fx - * - % - 


Differentiation yields the following results. 


p %* ■ ft 1 + V ^ 




■ h T ['- ?! $i T > ^ 
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”vi; • ■’ 


or 






3P 


A"! = ' 


• r * T 

X) = [ P X > 

Z T ] 

(6-2.2) 

: the partials are 

produced in like manner but for 


: results will be 

printed 


- 



3<x r %> 

“3 % 


. - . 

• 

3 (x x , X x ) 

(6-2.3) 

3 (x,i) 

3 % 



■a Pi / • = f J- ( e . pe 

7^,2) «*. ^ ^ 


, T 


/ 3-1 3(x,2L> 


a (%>*!) 

3 s. 


(6-2.4) 


(6-2.5) 


For use in Equations 6-2.3 and 6-2.5 above, we need the partial derivatives 
of the instantaneous geocentric ecliptic state of the tracking station 
(x^ and x^) with respect to either its spherical or cylindrical coordinates. 
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The geographic position of a tracking station is defined in terms of 
radius (r), latitude (<J)) , and longitude (A) when spherical coordinates 
are being used and in terms of spin radius (r g ), longitude (A) and 
z-height (z) for cylindrical coordinates. For most Earth orbit applica- 
tions, the tracking stations are specified in terms of spherical coor- 
dinates . 

If G represents the sidereal hour angle of Greenwich at launch, t 
the universal time (U.T) at the launch epoch, t the current U.T. after 
the launch epoch, and co the sidereal rotation rate of the Earth, then 
we have the geocentric ecliptic state of a tracking station being 
given by 


E 


for position, and 


I r cos <{> cos (A + G + oj (t-t Q )) 

1 r cos <J> sin (A + G + u* (t-t Q )) 

r sin 


*1 


E 


-tor cos <(> sin (A + G + <*> (t-t Q )) 

+cor cos <J> cos (A + G + to (t-t )) 

o 

o 


for velocity. If the geographic station location is specified in terms 
of cylindrical coordinates, then the geocentric ecliptic state is given 


by 
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1 


r cos (fl+ G +«*> (t-t )) 
s - o 


*1 


r sin(fl+G+co(t-t)) 
s o 


and 


-to r sin ( A + G + «*> (t-t )) 
s o 


*1 


+ 60 r cos (A + G + 6£ (t-t )) 
s o 


0 


The matrix E, appearing in the above expressions, is the 3x3 trans- 
formation from the geocentric equatorial to the ecliptic Cartesian 
frame. The elements of E are assumed to be constants even though 
they depend on the obliquity of the Earth's equatorial plane to the 
ecliptic plane. In other words, the temporal variation of the 
obliquity is assumed to be negligible over the duration of missions 
for which the program is used. 

First looking at the partials of and x^ with respect to the 
spherical coordinates, = (r, A ), we have 

cos cos 0 

cos sin 0 

sin <)) 

and 



-r sin § cos 0 -r cos <J> sin 0 
-r sin 0 sin 0 r cos cos 0 
r .cos $ 0 

-l 


0 
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- cO cos (J) sin 0 tor sin (}) sin 0 - cO r cos ^ cos 0 

cO cos <|> cos 0 Or sin <J) cos 0 - co r cos 4) sin 0 


where 0 is defined as 


0 = * + G + co (t-t Q ). 


Next we give the partials corresponding to those written above, 
but for the cylindrical coordinates, = (r g , X > z )« They are given 


cos 0 


sin 0 


-r sin 0 
s 

r cos 0 
s 


- go sin 0 


cos 0 


- to r cos 0 

s 

- tJ r sin 0 

s 
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Three-way range and range-rate are measured with one station on the 
DSN uplink and another station on the downlink. Three-way data may 
be processed by itself, simultaneously with conventional two-way data, 
or as differenced two-way minus three-way data, also known' as QVLBI 
(quasi-very long baseline interferometry). Three-way data types are 
modeled as the sum of the two-way types plus a timing error term for 
ranging and a frequency bias term for range-rate. 

<*3 ■ *1 + P 2 + cAt 

» • • A r 

^3 = + ^2 + C ~f 

where At is the timing error, c the speed of light, and ^/f the 
frequency bias term which results from drift error between the ‘frequency 
standards at the two separate tracking stations. The sensitivity 
partials for the three-way data types are formed by adding the par- 
tials computed for each station individually. The c At and 
c A ^/f terms are treated either as biases or part of the white 
noise term. The differenced data types are modeled: 

AP = ( 5 1 - e 2 - cAt 

a£ = - e 2 - c Af /f 

The partials for the differenced data types are formed by differ- 
encing the individual partials, with the following exception. Since 


j 
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SAe 

c)x 



A 

and ^ and 
Earth orbital 


A 

P 2 are very nearly equal (as are and for 

missions, we use the following substitutions 


^ x = x 2 " Aj 

I A A 7 

?1 + ? 2 1 

A T 5 ~ 

1+ Ci e 2 

Cj - f 2 ” [a* - aP P 2 ]/p! 

2&£. - l*t - [A* -APP 2 l/» (6-2.8) 

3 X 3 i 

For Earth orbital missions, one of the key navigation instruments 
is the infra-red horizon sensor (H.S,). The basic H.S. measurement 
consists of four angles which are meausred by optically scanning the 
figure of the Earth in two orthogonal directions. Two angles are 
evaluated in each direction and measure the apparent separation 
between a reference measurement axis and the optically sensed contracts 
with the horizon. Since these data determine the angular diameter and 
direction cosines to the Earth, a single H,S, observation is sufficient, 
in theory, to fix the S/C position. The only additional requirement is 
an accurate estimate of the S/C inertial attitude. 

For the purpose of mathematically modeling H.S. observations, the 
set of four angles is reduced to three equivalent angular measurements. 
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©^ © 2> and © 3 . These new angles define the position vector, £, to 
the geocenter as seen in the measurement coordinate system. (See 
Figure 6-2.) ©^ is measured in the e^, e^ plane and is the angle 

between the ^ axis and the projection of £ into that plane. © 3 is 
analogous to ©^ except it is measured in the e^, plane, relative 

to the e„ axis. ©„ is just the apparent angular radius subtended by 

J ti 

the Earth. Thus, the observation equations are given by 



Figure 6-2. 


Horizon Sensor Measurement System 
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Since the attitude of the vehicle is nominally assumed to be known, the 
orientation of the e^, e^ 3 axes is defined such that points - 

vertically down, is directed along the local horizontal in the direc- 
tion of motion, and forms a right-handed orthonormal traid. With this 
definition for the e^, system, an arbitrary vector can ^be con- 

veniently mapped into an inertial representation by multiplying the 
vector by the transformation matrix, R, defined by 



I 

I 

» 

1 

\ 


v x r_ 

i- x -i 


i 

i 

i 

\ 



and £ and v are the position and velocity vectors evaluated on the 
reference trajectory at the observation point. 

In GODSEP, the quantities essential to the error analysis are the 
observation partials which relate the measurement residuals to the geo- 
centric state of the S/C. To obtain these partials, the observation 
equations are differentiated with respect to £. This gives 




H Sj£ 
P 


where 
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To relate these measurement residuals to 
rotation matrix discussed above, i.e., 

=-H , 

P 

where & r is a position deviation relative 
negative sign appears here because ]> is in 



the geocentric state, we use the 


to the reference position* The 
a direction opposite to r_. 
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The following definitions are used in the azimuth and elevation 
angle partials (See Figure 6-3). 


at = S/C azimuth, measured positive from north toward east 


3 

e 


A 

e 

?, 


X 

X 


S/C elevation 

S/C range vector from station 

Unit vector in £ direction 

Projection of P onto plane normal to x g 

Geocentric equatorial S/C position 
Heliocentric ecliptic S/C position 


x 


A 


X 

A 

W 


s 


A 

u 


E 


Geocentric equatorial station position 

Unit vector in x direction 
— s 

Unit vector orthogonal to x and pole (local east 

s 

from station) 

Unit vector orthogonal to and w (local north from 
station) 

Transformation from equatorial to ecliptic coordinates* 


For simplicity, all azimuth and elevation partials are derived in 
geocentric equatorial cartesian coordinates and then transformed to 


ecliptic. 
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sin $ 


* T Z 

• \ e 


11 , i 

3 X Q cos |J 


[ A A X 

x s - sin(J e ] 


(6-2.9) 


_ 1 

3x s |Ss| cos P 


[ § - sinfl x ] - 


M. 

c> X 


(6-2.10) 


,(S /^ - */»*. 3 -/ 3£ 


(6-2.11) 


Again referring to Figure 6-2, the projection of onto tf will have 

gnitude cos ^ sin at , or 


.ma 


sino(. = sec £ £ p 


If ■ can * can P + 


3 x. 


a, = | secc^ sec w - tanat £ i /e 


( 6 - 2 . 12 ) 


= tan<*.tan£ ig + 


where 


b = 


” [sec*/ sec Q (w x C y "« y f x )]w / (6-2.13) 


^ A ^ ^ 

fx, f y, W^, are components are w and 
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c>*/ 



For use in Equations 6-2.11 and 6-2.13 above 


(6-2.14) 



d * s 

9 (r s »l,z) 


X, /? -x* 0 

1 / s 2 

s s 

x 2 / X]L o 

S S s 

0 0 1. 


(6-2.15) 


where x- and x„ are components of x . Finally, the partials 

Jl / S 

s s 

computed in equatorial coordinates must be transformed into ecliptic 


>(«.*) = ISJL ft ) . lA. 

*a *c 92^ 


^ 9 ) = 9(°4> £ ) . E T (6-12.16) 

92^ t 2L 


For vehicle based optical measurements other than the horizon 
sensor observations, we use the following defintions: 

£ = spacecraft geocentric Cartesian position, 

A 

= unit vector of star direction cosines, 

R = radius of the Earth, 
e , 

f 

V* = star-Earth horizon angle, 

£ = apparent planet diameter measurement - angle substended 
by planet disc at the spacecraft, 
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z = 3x1 null vector. 


Star-planet angle partials: 


c 


/ 


ns 


cos 
^ r 


* = 


a ? 


f sLy W - <=°sYf > 


Apparent planet diameter partials: 


sin 1/2 £ = ~ 

r 

r 

= -4 tan I/ 2 S P 
fcf I I 

For any data type which has a bias, 


y - H % X + b 


3y 


/ = 1.0 

' * b 


(6-12.17) 


6.4 FILTER 

After the knowledge covariance has been propagated to a measurement 
time, and the observation matrices are computed, the 0D filter can per- 
form its function of estimating the set of solve-for parameters (non- 
deterministically) and updating the knowledge covariance accordingly. 
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There are two types of filters available, Kalman-Schmidt (K-S) and 
weighted least squares (WLS), plus capability for a third filter, to be 
established by the user. K-S is the most commonly used filter because it 
treats consider parameters in a realistic fashion. 

The filter updating process requires computation of several matrices. 
First, the propagated estimation error (knowledge) covariance P at the 
measurement event time is 
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P = 


xs 


X 

: T P 
xs s 


u 


xw 


xw 


P c 

V vw 
T 

C 1 P 
vw w 


where P , • P are the covariances of the S/C state, ignore para- 

x w 

meters, respectively, and C^, •••, are the cross-covariances between 

appropriate augmented parameters. The observation matrix H defined in the 
previous section is 

* H 


H = 


H 

s 

H 

u 

R 

v 

. H i 

l w J 


The measurement residual matrix J is defined as 

J = HPH T + R 

where R is a diagonal matrix containing variances of the measurement white 
noise. For example, a simultaneous 2-way/3-way range measurement would 
look like 

1 ^ 

'2R 

0 


R = 


<*2R + « 2 


>3R J 


where & 2 is the 2-way range noise variance and is the additional 

2R 

3-way range noise variance due to timing synchronization. For a single 
star-planet angle measurement, R would be a scaler 

r - <r 2 + o.V 
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2 _2 

where and 0^ are the optical resolution and planet or body 

center finding noise variances, respectively, and r is the spacecraft 
range to the planet. 

The updated covariance (P + ) after the measurement has been 
processed (and in theory after the state estimate has been updated) 
is in general 

P + = (I - KH) P (I - KH) T + KRK T (6-3) 

where K is the filter gain. 

) 

KALMAN-SCHHIDT 

The filter gain for K-S is straightforward 
K = PH T J -1 

Since only estimated parameters can be updated by the OD process, 
the entries of K corresponding to consider and ignore parameters 
are zeroed out, that is, K = [ \ K g 0 0 0 ] T . The updated 

covariance is then formed by Equation 6-3. 


WEIGHTED LEAST SQUARES 

The sequential, or recursive weighted least squares (WLS) algorithm 
implemented in GODSEP is equivalent to a batch WLS filter i^f there is no 
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process noise* Since process noise is a significant part of low thrust 
analysis, the WLS filter must be used recursively, because it has no 
batch equivalent. The sequential WLS consider filter acknowledges con- 
sider parameters only in the covariance update (Eqn. 6-3) and not in the 
gain matrix calculation. Therefore a set of reference covariances 
for the state and solve-for parameters must be maintained at all times. 

This set also represents the filter analysis as it would be in non-consider 


form. 

Thus, the WLS filter computation requires three operations: (1) pro- 

pagation of the reference covariance to the measurement event, (2) com- 
putation of the filter gain and (3) updating both the reference and 
knowledge covariances. 

(1) The reference covariance (P) consists of 


A 

P * 


xs 


xs s 

and is initialized at the a-priori values. Thereafter, it is 
propagated from one measurement to the next by 

p = 4 p 
*k+l 9 k p 

where i is the augmented transition matrix corresponding to the 
x and s parameters. P is computed in parallel with the actual 
knowledge covariance P. 

Given P at the measurement event, the WLS filter gain is 


a a T* * *T -1 

K = P H (H P H + R) 


where H - 


m 

H 
_ H 


( 2 ) 
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(3) The reference covariance is updated, after measurement pro- 
cessing, by 

$ + = (I-KH) P 

and the knowledge covariance P is updated by Eqn. 6-3. 

Thus, at measurement events, the 0D filter updates the knowledge 
covariance to simulate taking a tracking measurement, processing the measure- 
ment in an orbit determination algorithm, estimating desired parameters and 
reducing (or updating) the knowledge uncertainty to reflect this new infor- 
mation about the trajectory. 

6.5 GENERALIZED COVARIANCE 

The function of any filtering algorithm is to process available measure- 
ment information and produce a best estimate of the spacecraft state and 
any parameters that are being solved-for. Best is usually defined in a sta- 
tistical sense, such as the minimum variance processes used in differing forms 
in the weighted least squares and Kalman- Schmidt filters. But in practice, 
filter performance is dependent on how well the assumptions used in the 
filter definition approximate real-world processes, because all error sources 
cannot be modeled, nor can those that are modeled ever be modeled exactly. 
Therefore, each filter must be evaluated not only on its ability to produce 
small error covariances in the resulting estimated state, but also be as 
insensitive as possible to errors in its model assumptions. 

Generalized covariance error analysis is a useful tool for studying 
filter sensitivity. For generalized covariance studies, two sets of know- 
ledge errors are carried during the orbit-determination process. Assumed 
knowledge uncertainties are those generated by the filtering algorithm 
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according to the mathematical model and all the assumed errors input to 

„ l 

it. True knowledge uncertainties represent the effect the filtering 
algorithm has on actual state estimation when the real-world error sources 
are not the same as those assumed by the filter. Evaluating filter sensi- 
tivity to a model assumption involves comparing the resultant effect on 
assumed and true uncertainties of a modeling mismatch between the filter 
and real-world uncertainties. This modeling mismatch is accomplished in 
one of two ways; true a priori uncertainties may be set at levels other 
than assumed levels or the true state may be augmented by a vector of ignore 
parameters — parameters whose uncertainties are recognized by the true 
covariance analysis, but which are completely ignored by the assumed filter 
analysis. 

The filter that is least sensitive to a model mismatch is determined 
on the basis of two criteria. First, which filter yields the smallest true 
estimation errors. Second, for which filter are the true errors most closely 
approximated by the errors predicted by the filter covariance analysis. Thus, 
given a mission with a specific set of model mismatches, if two different 
filters produce equivalent true errors, then the superior filter is the 
one whose assumed errors are closest to the true ones. Similarly, if the 
resultant assumed errors of the two filters are equivalent, the superior 
filter is the one with the least true estimation uncertainty. Generally, 
qualitative judgments are required because several sets of mismatches must 
be studied, and the relative performance of the filters may vary. 

In error analysis, generalized covariance is a filter sensitivity 
study tool that is normally available only in a simulation program. It is 
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accomplished in G0DSEP with a minor increase in core and computational 
time compared to a full simulation and has the additional advantage of 
generating ensemble true state statistics rather than a single sample as 
in a simulation. The only disadvantage of generalized covariance is that 
it uses the same linearized dynamic and observation models as the assumed 
filter analysis, and can therefore not study problems that arise from 
nonlinearities. 

The actual operation of generalized covariance in G0DSEP requires that 
a standard error analysis be run first. The filter gains, associated with 
the assumed knowledge uncertainties, are stored on disc or tape. Now the 
error analysis with all the same measurement events is repeated. Only this 
time, a-priori uncertainty levels and measurement noise are modified, and 
ignore parameters are added, to the extent of desired mismodeling. At 
each measurement event, Eqn. 6-3 is applied to what is now the true know- 
ledge covariance using the appropriate stored filter gain. The true 
covariance analysis thus proceeds in analogous fashion to the previous 
assumed covariance analysis. Obviously, many mismodeling conditions can 
be studied with the same filter by repeating the generalized covariance 
analysis. 

6.6 guidance 

Although the knowledge covariance is modified by measurement events, 
the control covariance, which represents the ensemble of actual deviations 
from the desired or reference trajectory, will grow without bound. The 
only process which will reduce the control covariance is a guidance event 
that is, the design and execution of a trajectory correction, either 
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impulsive AV or low thrust. 

Low thrust guidance represents an update of the nominal thrust con- 
trols (magnitude, direction and cut-off time)* In terms of system cost 
and efficiency, it is better to use the existing low thrust propulsion 
system for guidance than to add auxiliary means, for example high thrust 
chemical engines to produce impulsive iv. Of course, certain problems 
inherent in low thrust propulsion, in particular terminal controllability, 
may force the addition and use of an auxiliary chemical propulsion system. 

In mathematical terms, given a trajectory state deviation, £Xo = 

$ x (t Q )j where t* is the guidance epoch, we wish to null out the effects 
of S x by making a bias type correction $ u to the nominal thrust controls. 
To be efficient, the correction is applied over some finite interval 

r t , t I such that the target error £t , caused by $ x , at some final 
time (t f ) is removed. For linear analysis, we seek the guidance matrix P , 
such that 

8n - r 


The linear ensemble of thrust control corrections is then 


U = E [ S il £ j± T J 

u = r Efs^ sx/] r 


(6-4) 


In G0DSEP, the trajectory error ensemble E £ £ x q Sx^jis the control 

covariance (t Q ) . Using represents a pessimistic sizing of the 

thrust corrections because only the known trajectory error (not to be 

confused with the knowledge error) can be removed. The known error 

generally corresponds to the control error as long as P (t )>> P (t ) 

co k o 

where is the knowledge covariance. 
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To compute the guidance matrix f we first compute the sensitivity 


matrices 


or. 


S 


V 


s (t f ,t o ) 


3 T (t f ) 


31 (t f ) 
^ x (t o ) 


rjaicv j 

l J 

l J 


♦ <t £ ,t o > e (c c ,t 0 ) 


where the first matrix in S and V is formed by numerical differencing 
and the second two matrices ((|> and 0) are obtained from transition matrices 
generated by the trajectory propagation routine (Section 6.2). If variable 
time of arrival is desired, the control array Ju is augmented with the 
arrival time and the S matrix is augmented by x (t f ) , relative to 
the target. 

The guidance matrix can now be defined by 

P = -w 2 s T [ SW 2 S T 1 V If N > N (6-5.1) 

u l u J u — T 

P = -[s T W*s] 1 S T W t 2 V If N t > N u (6-5.2) 


where and are the number of parameters in the control and target set, 
respectively, and and are diagonal weighting matrices for the control 
and target parameters, respectively. The first form of P, 6-5*1, reflects 
a minimum quadratic control correction and the second, 6-5.2, corresponds 
to minimum quadratic target error. 
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Generally, there are more control parameters than targets with the 
exception of two cases: (1) terminal approach to the target where con- 

trollability drops off rapidly; and (2) the application of control 
constraints which effectively reduces the set of available controls. 

If constraints on the control corrections are imposed, then P must 
be suitably modified. First, the unconstrained P is computed along with 
the ensemble unconstrained control corrections, 


u = r V r T (6-4) a s ain 

2 

Each diagonal component of U, ffu , is compared against its constraint 
value ( Su^) 2 . If (f u is greater than S u^, then the appropriate 
row of T is scaled by The total control set (and guidance 

matrix) is then separated into two subsets: unconstrained controls, £ U| 

and P , and constrained controls, £ u^ and P^. 



The new control corrections are computed with (6-4), (actually only the 

remaining unconstrained controls are computed ) the test for constraints 

are made, P is modified again, and the entire process is repeated until 

all constraints are met, or there are no more controls left. The guidance 

corrections are executed (figuratively) at time that is, are uplinked 

to the spacecraft, but apply over the entire guidance interval [ t Q , . 

Execution of the guidance updates causes the control covariance to 

diminish from t to t whereupon it begins to grow again. Guidance accuracy 
o c 

is measured by how much the control covariance can be reduced at t which 
depends upon how well the thrust corrections were designed. The limiting 
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accuracy of maneuver design is the knowledge error or covariance at t Q . 

Thus, the post maneuver control covariance at is the propagated 
knowledge covariance (as per Eqn. 6-1) from t Q to t c> 

F+(t c > = ♦(t c ,t 0 ) Ct 0 ) <t c ,t 0 ) +Q(t c ,t 0 ) (6-6) 

In G0DSEP, to denote guidance execution, P+(t Q ) is set equal to 

P (t ). This is equivalent in effect at t to applying Eqn. 6-6. 
k o c 

However, it means the value of P* in the guidance interval is not valid. 

This is a relatively minor problem compared to the reduced burden on 
computational storage and logic. 

One exception to setting P* = P k occurs when there are more 
targets than available controls, which often happens when control 
constraints have been activated. In this case, there will be some non- 
zero target error that was not removed by the guidance corrections. This 
implies that not all of P c was removed. Hence, the post maneuver control 
covariance must include the residual state error, 

P+ = P k + ■[ v T [w T ] [ V + s r] | P c " | v T [ w T ] [ v + s flj T 

As long as no more guidance events are executed between [t Q ,t c } , 
updating the control covariance at t^ is theorettically no different 
than using Eqn, 6-6 at t • However, if another guidance event is scheduled 
in the burn interval, say at t^ then a somewhat different logic is applied 
Co size the correction, Xt is assumed that the first guidance event between 
t and t is a primary maneuver. Subsequent guidance events in this inter- 
val are considered to be vernier maneuvers, representing refinements of the 
thrust corrections computed in the primary maneuvers. 
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For vernier maneuvers , 0 and (f are redefined for the vernier burn 
interval and P is computed as in Eqn. 6-5. The guidance updates are 
computed using Eqn, 6-4 with the knowledge gained since the primary 
(or previous vernier), that is, 

E [ *£[]- P X> - W 

Recall from the previous discussion that P^(t^) is usually the propagated 
knowledge from the previous guidance event. 

A measure of guidance effectiveness is the estimated target error 
before and after the maneuvers. 

before guidance correction, E [$T S T T J = V P (t ) V T 

f rjy -j rri 

St St } = vp(t)v. 

~ 4 CO 

This simple measure assumes, of course, that no further dynamic error 

will occur from t to the target time t £ . 

o f 
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An important part of the guidance and navigation process is 
the time interval from the last navigation measurement used for 
guidance design to the actual time of guidance implementation 
(t Q ) . The time interval (or delay) is necessary for ground pro- 
cessing of all previous measurements, estimation of the S/C state, 
designing the thrust updates to correct the trajectory, and 
execution of the updates. Typical intervals are 3 to 15 hours, 
and are usually critical only in the terminal mission phase where 
trajectory controllability (with respect to thrust controls) 
diminishes rapidly. This time delay is user specified for each 
guidance event. 


Impulsive &V guidance is very similar to low thrust except 

for a zero burn interval (t = t ). The delta-velocity is treated 

o c 

as if it were a control correction § u, that is, 


Ai = P 

To compute P * the sensitivity matrix S is first partitioned 

into position and velocity submatrices, 

« 

s = [ a : b ] 


*>=0 


31(t f ) 


where A 


and B 
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If the target vector T has only two components, e.g., B*T and B*R, 
then 1 AY l is minimized and 

P - [ -b t (bb t ) A [ -b t (bb t ) B ] 

If T has three components, e.g*, the position vector at t^, then 

r - [ -b- i a i -i ] 

The computation of ensemble velocity correction, U in Equation 6-4, 

follows directly. As in low thrust guidance, the pre-maneuver control 

covariance P (t ) is used to size U. 
c o 

u = e[&vay t ] = r p; (t 0 )p T 

Execution errors related to low thrust control updates are neg- 
lected because they are second order effects compared to thrust error 
associated with the nominal thrust profile. However, AV execution 
errors are taken into account because impulsive maneuvers often occur 
during ballistic or coasting portions of the mission and can represent 
a significant contribution to trajectory error. In order to compute 
&V execution errors, the most probable &V is first determined by 
the Hof fman- Young approximation (Reference 8 ). Let A ^3 

be the eigenvalues of the AV covariance, U, and oi be the largest 
eigenvector of U, define 

A = > 1 +X 2 + A 3 

B=XiX 2 + ^1^3 + ^2^3 




8? 


then the probable AV is 


e[av] 



AVi 

av 2 

av 3 


Now the 3 x 3 AV execution error covariance Q is composed of 

*q - AV , 2 <r 2 + AV, 2 P 2 K 2 + ^i 2 av/ V (6 

11 1 m p 2 

~xy 

Q 12 -Q 21 - av, AvJrr 2 - P 2 O ', 2 + AV , 2 <T, 2 ] 

1 ” ^7 J 









2 
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2 

(T = out of ecliptic (z) pointing variance. 

v 


As in low thrust guidance, the post-maneuver control covariance P c 
(t ) is set equal to the knowledge covariance P, (t ) corrupted by 

O K O 

the execution errors. 




( t o ) + 



Pre and post maneuver target uncertainties are computed in equivalent 
fashion to low thrust guidance. 
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7.0 TRAJECTORY SIMULATION - SIMSEP 

The trajectory simulation mode SIMSEP has been designed to 

deterministic analysis of ballistic and low thrust missions* 
Computationally, SIMSEP imitates "real” trajectories in the presence 
of a wide variety of environmental and system uncertainties. A pri- 
mary objective is to deduce expected or probabalistic behavior of 
the real mission by studying a relatively small subset of simulated 
missions. 

The purpose of this section is to discuss the key analytic 
concepts in SIMSEP. This will be done in two parts; 1) by discuss- 
ing the principal algorithms, and 2) by outlining the basic compu- 
tational structure. Although many algorithms used in SIMSEP are 
similar in function to algorithms used in T0PSEP and G0DSEP, their 
specific applications here warrant an extended discussion of their 
underlying theories. 

7.1 Program Scope and Methods 

Before proceeding with a step-by-step description of the 
algorithms and computational structure, it is worthwhile to compare 
the essential similarities and differences between G0DSEP and SIMSEP. 
Unlike the error analysis mode which works exclusively with error 
ensembles and a reference trajectory, the simulation mode actually 
formulates many discrete examples of the "real world" or "actual" 
trajectory. Each of these is propagated, in a deterministic sense, 
by the same trajectory integrator, integrating the same equations of 
motion. However, many variables and parameters appearing in these 
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equations are subjected to random alterations, corresponding to 
discrete uncertainties. Hence, each deterministic simulation of a 
mission is different according to .the effects of the sampled errors. 

Operationally, SIMSEP does not sample each error in succession 
and propagate trajectories with just one error source active 
time. Rather, all are initialized by differing amounts for each 
actual trajectory. Thus the averaged effect of all error sources 
acting in concert can be estimated by repeating the mission simula- 
tion process a sufficient number of times.. This is the essence of 
the Monte Carlo method and is the basis of the simulation approach 
to determining how trajectory nonlinearities and uncertainties can 

affect the G&N process. 

Perhaps the beet example to illustrate the fundamental differ- 
enees between the methods used In G0DSEP and SIMSEP Is the problem 
of propagating, or mapping, an error eovariance from one point to 
another along the referenoe trajectory. It will be recalled that 
the principal method for propagating a covariance in G0DSEP is by 
the state transition matrix mapping, namely. 




lk+1, k P k *k+l, k 


where $ k+1 k represents the state transition matrix and P k and 
P k+1 are covarianceS at h a " d 'fchl* respectively. When there is 
dynamic process noise, an effective process noise matrix, Q k+lk . 
is also added, (See Section 6.2). The state, transition matrix is 
generally computed simultaneously with the trajectory by integrating 
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the variational equations. However, the variational equations are 
based on linearized expansions of the differential equations of 
motion and neglect all second and higher order terms. Hence, a 
state transition matrix mapping of covariances must theoretically 
be limited to mapping covariances within the envelope of linearity 
surrounding the reference -trajectory. Rarely is this assumption 
true at all points along an interplanetary trajectory, especially a 
low thrust trajectory. Covariances propagated by this means are 
subject to error whenever a region of significant trajectory non- 
linearities is encountered. 

On the other hand, the method for propagating a control 
covariance in SIMSEP is not plagued by these effects, although it 
has its own peculiar shortcomings. The SIMSEP approach to this 
problem relies on the Monte Carlo method where a multitude of 
sample trajectories are propagated between the two time points in 
question. The trajectory state vector data at are processed 

and accumulated in such a way that the covariance can be reconstructed 
by standard statistical calculations* Hence, SIMSEP maps a covariance 
as a statistical ensemble calculated from many data points and not as 
a simple mathematical entity like G0DSEP* 

Therein lies the primary drawback to the Monte Carlo method and 
to the use of SIMSEP for general G&N analysis. Although the Monte 
Carlo method will, in theory, converge to the exact covariance, the 
rate of convergence tends to be extremely slow. For accurate 
statistics, inordinately large amounts of computer time are often 
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necessary to perform the many trajectory propagations. 

Although both G0DSEP and SIMSEP are preflight mission analysis 
programs used to identify the general G&N subsystem characteristics, 
they are usually used at differing phases in the development of an 
overall systems analysis. For the purposes of a preliminary systems 
design, a G0DSEP analysis is the most cost effective means of evaluat- 
ing the basic G&N subsystem requirements. As such, SIMSEP is generally 
relegated to verifying the linear analysis results, and only in the 
advent of serious nonlinearities is the simulation mode called upon 


for more extensive studies. 

7.2 Definitions and Concepts 

The first important concept in SIMSEP and common to all MAPSEP 
modes is the reference trajectory, denoted by X^ = X^t). The 

reference trajectory is computed under some set of "reference 
integrating conditions" to satisfy desired targets at mission's end. 
Moreover, represents a deterministic solution to the equations of 
motion for the assumed dynamic and systems models. For SIMSEP, the 
initial state and reference integrating conditions, i.e., ephemeris 
parameters, thruster characteristics, etc., are read as input since 
it is assumed that they have already been computed as output from a 
T0PSEP analysis. 

A second quantity important in SIMSEP and common to G0DSEP is 
the control error covariance, P » Generally, an a priori control 
covariance is defined at injection (or at the starting point of the 
mission being studied). This matrix mathematically describes the 
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distribution of real state errors relative to the initial reference 

trajectory state. In SIMSEP, it is implicitly assumed that the 

probability distribution of these errors is Gaussian with zero mean. 

Once an a priori control has been given, it is randomly sampled to 

form an error vector, 8 X. , which corresponds to a deviation of the 

A 

actual trajectory state relative to the reference. At the same 
time, error sources associated with the host of other dynamical 
and systems uncertainties are also sampled to create the so-called 
M real world integrating conditions". For the actual trajectory 
state vector, this procedure may be written as 

h ‘ h* s h 

where 8 X* is a deviation obtained by sampling P . Utilizing 
these integrating conditions, the actual trajectory state, X^, is 
propagated from point to point as a discrete example of an actual 
trajectory. 

The third- v cri tical variable used in both GODSEP and SIMSEP is 
the knowledge error covariance, P R . This matrix is propagated in 
the error analysis from measurement to measurement where it is 
systematically updated according to the filtering algorithm. In 
SIMSEP, instantaneous evaluations of P^ are input at each guidance 
event and are left unchanged throughout a given run since there is 
no explicit orbit determination process modeled. However, the 
knowledge covariance, like the control covariance, is sampled to 
formulate an error vector. This error vector determines the error 
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in the estimated state relative to the actual trajectory state, and 
is used to compute an estimated state vector by 

h = ^ + ^ 

where e is the sampled error vector from P^. If other parameters 
are estimated during the orbit determination calculations, they are 
included as augmentation parameters. These too are sampled in order 
to formulate a set of "estimated world integrating conditions". 

With each of these key quantities having' been defined, it is 
worthwhile to mention why each is important in a simulation run. 

The reference trajectory, for example, serves to define the mean 
for all actual trajectories, as well as defining the reference target 
conditions used during guidance. The actual trajectory is, of course, 
the mathematical representation of the real motion and is carried 
from event to event until the final target is reached. On the 
other hand, the estimated trajectory is used exclusively for re- 
targeting the actual trajectory back to the desired targets and is 
computed only during guidance. 

7.3 Guidance 

One of the principal purposes of SIMSEP is the detailed exami- 
nation of nonlinear trajectory effects, especially as they bear 
upon the guidance problem. In this section, the fundamental concepts 
underlying linear and nonlinear guidance will be presented, and the 
implementation of these concepts into algorithms will be discussed. 
Beforehand, a careful distinction between targeting and guidance 
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must be drawn since MAPSEP has algorithms for both and since many 
of the basic steps and operations are similar. Whereas a targeting 
problem is solved by formulating an entire control strategy for a 
complete mission, a guidance problem assumes that a solution to the 
targeting problem has already been found. Furthermore, the current 
trajectory which is to be corrected is assumed to be in a "close 
neighborhood" to the original reference solution. Hence, the con- 
trol changes computed by a guidance law are expected to be small 
refinements to the original controls, even in the presence of non- 
linearities. 

7.3.1 Linear Guidance 

The linear guidance option in SIMSEP is analogous to the guid- 
ance used in G0DSEP except that it applies to a discrete trajectory 
error as opposed to an ensemble of errors. For both modes, the 
linear guidance matrix is the same. To compute a guidance matrix, 
a sensitivity matrix is evaluated between the point of the guidance 
event and the target about the reference trajectory. This matrix 
of linear partials relates control changes to target deviations and 
is used to map estimated trajectory errors, SXg, into control 
updates, according to the guidance laws: 

& V - T S X^. , for impulsive corrections, and 
2) _Au = f _Sx E , for low thrust. 

In spite of its overall simplicity and ease of implementation. 
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the advantages of linear guidance are off-set somewhat by its draw- 

t 

backs. First, the trajectory error, S X^ , mast lie within the 
envelope of linearity for this to be a valid method. Whenever this 
is violated, the resulting guidance correction can be invalidated. 
Furthermore, a linear guidance correction is executed without itera- 
tions. In fact, with this guidance there is no direct evaluation of 
a control correction's effectiveness in reducing target error. Only 
if the updated trajectory is propagated to the target can the result 
ing target error be determined, and even then there is no recourse 
for making further corrections if the original correction is 
ineffective . 

7.3.2 Linear Impulsive Guidance 

The essence of impulsive guidance is founded on the mapping 
relations which propagate arbitrary linear deviations relative to 
some known trajectory into new deviations at some later time. 
Clearly, this is a property of the state transition matrix which 
maps a six component state deviation, evaluated at into a 

new deviation, S ^ J r at t k+1> by the equation, 

Men ' W-^4 


If t is the target time and t, is the time of the guidance event, 
k+1 k 

then $ can also map state vector changes, like an impulsive 

Ik+l,k 

velocity correction, into state changes at the target. 

However, in most analysis the actual target conditions are 
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specified in terms of target variables such as geographic coordinates, 
Keplerian elements, etc., instead of X, Y,...., Z state coordinates. 
Fortunately, the target variables are functions of the final trajec- 
tory state and it is possible to generate a differential transforma- 
tion of the form 


lx - ^ 14+1 <7 - 2) 

which transforms differential coordinate changes into target variable 
variations* In the above equation, represents a matrix of linear 
partials of the form 


'V\ = ^ * T 1» T 2» ’V 

\ " * (X, Y, Z, X, Y, Z) 


k+1 


where there are n-target variables, T^, T ... 
components. By substituting Eq. 7-1 into 7-2, 
state changes at into target changes at Vi 


It - 


•k+l , k 


Sx, 


(7-3) 

.. T and six state 
n 

a relation for mapping 
is obtained, that is. 


Performing 

~\ $H-l,k 


the indicated matrix multiplication and 
with N, the equation becomes 


replacing 


&T = N _SX k 


(7-4) 


where N has dimensions (n x 6) . 

With this background, the impulsive guidance problem can be 
stated most simply as a determination of a velocity change, A V. 
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which, when added to li, nulls the target error, ST . Again 
writing Eq. 7-4 but in a partitioned format, 

St = N(l) |r k + N(2) Sy k , 

it is recognized that h T can be made zero by adding some appropri 
ate AV to Sy_ k > i.e. 

N(l) % r fc + N(2) (&v k + &T) = 0. 


For the case of three-variable impulsive guidance, i.e., three 
unique targets, the solution for AV is given as 


4V = - N(2) _1 K(l) 


provided N(2) is nonsingular. Note that this can be re-written as 


= V-NG)' 1 n(1) ; - 


(7-5a) 


to 


or 


4y - v 


(7-5b) 


the desired guidance law. 

For the case where there are two target variables instead of 
three, the problem has more controls (3-velocity components) than 
end conditions and a generalized inverse which minimizes the magni 
tude of the velocity correction is used according to 
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&L - - N(2) T [n( 2) N(2) T 1 N(1) — k 


- N(2) T ^N(2) N(2) T *| 1 iv^ 


where N(l) and N(2) are non-square matrices with dimensions (nx3). 
Again this relation can be re-written as 



(7-6a) 


(7-6b) 


Algorithms based on Eqs. 7-5a and 7 -6a are the basis of the 
linear impulsive guidance contained in subroutine LGUID. The guid- 
ance matrix, T"* , for either the two or three variable cases, are 
computed as outlined above and the state vector deviation, £ X , 

K 

is calculated as the error in the estimated trajectory state relative 
to the reference, namely, 

- SX E * h - V 

evaluated at the guidance event. 

7.3.3 Low Thrust Linear Guidance 

The low thrust linear guidance law has the same format as the 
impulsive law except control changes are made to the vector of low 
thrust control variables, Another difference is that the low 

\ 
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thrust acceleration acts slowly to bring about state changes; hence, 
the low thrust linear guidance matrix operates in an integrated 
fashion over a fixed trajectory segment to redirect the motion. 
Otherwise, low thrust guidance is simply an extension of the basic 
methods discussed above. 

The most complex part of computing low thrust corrections is 
the determination of the guidance matrix, T"* ♦ This matrix depends 

not only on the trajectory dynamics between the maneuver point and 
the target, but it also depends on the trajectory response to con- 
trol changes, i.e., controllability. As before, the first step is 
to integrate the reference trajectory from the guidance point to 
the target, evaluating the augmented state transition matrix. In 
SIMSEP, the transition matrix is computed by integrating the vari- 
ational equations as was discussed in* Section 4-5. By selectively 
partitioning the transition matrix, the requisite sensitivity matrices 
relating state and control variable deviations to future state devia- 
tions are obtained. 

In terms of partitions in the augmented state transition matrix, 
state deviations at the target time, t^ +1 , are S* ven by 

- W + © <7 ' 7) 

where ® , , is a state transition matrix as defined in Eq. 7-1 

* k+1 , k 

and © is a matrix which maps control variable deviations into 
state changes at t^ + ^. S u in Eq. 7-7 corresponds to a set of 
thrust control biases. © can also be written as 
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0 = ^>(X, Y, Z, X, Y, Z ) m 

*U 2 f * • * » ^|» ) 


(7-7.5) 


and is seen to be (6xm) where m is the number of controls. Follow- 
ing the same line of reasoning as was given in Section 7.3.2, it is 
recognized that partials of target variable variations with respect 
to control variable changes are needed* Hence, Eq. 7-7 is multiplied 
by the transformation matrix (See Eq . 7-3) to obtain. 


St - 1 fk+i.k m n® ■£*• (? - 8) 


Therefore, the guidance problem is reduced to finding a An which 
when added to will make &T = 0. For convenience, it is assumed 

that & u is either zero or that it can be solved-for during the orbit 
determination; thus permitting Eq. 7-8 to be re-written as 




© Ajfe + 


^ i+i,k M 


0. (7-9) 


For the problem where the number of controls (m) equals the number 
of targets (n) and the matrix ^ has an inverse, the solution for 
A u in Eq. 7-9 is 


A n 



(7-10) 


where fe is the deviation of the estimated state vector relative 
to the reference at the guidance event. From Eq. 7-10 it is easy 
to see that the desired guidance matrix for this particular case 


is given as 
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- © _1 W 


and is a (6x6) matrix. 

The problem is complicated somewhat when the number of targets 

is less than the number of controls. In this case, a generalized, 

or pseudo-, inverse matrix operation is used, and the transformation 

matrix does not drop out. Nevertheless, a solution is obtained 

by determining the A*U- that makes S T = 0. Letting A = © and 

B = *>% in Eq. 7-9, a particular solution (out of the 

* kf 1 , k 

infinity of solutions) is given as 




i h- 


(7-11) 


This particular choice of AH. also minimizes the magnitude of the 
control change (See Section 5.3). Therefore, the desired guidance 
law can be written as 

Au = r sxg 

where T - - A T £aa T ] B. 

As before, the computational steps described here are imple- 
mented in LGUID. 

7.3.4 Nonlinear Guidance 

Nonlinear guidance parallels in many respects a targeting 
problem where an iterative, linear algorithm is used to determine 
control changes. The primary difference is that the trajectory is 
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assumed to be reasonably close to the reference. By reasonable, it 
is meant to imply that the algorithm should be able to redirect the 
s/c motion to the designated target by making a few iterations (three 
or four). In practice, a real trajectory can deviate widely from 
the reference and thereby require as many as eight to ten iterations 
before the guidance algorithm is able to compensate. For situations 
where convergence is not achieved after many iterations, the guidance 
is said to be divergent. The real criterion for qualifying a maneuver 
as divergent is somewhat subjective and established by the analyst. 

In some cases, an extremely slow rate of convergence on the part of 
the linear correction scheme can be attributed to non-adaptive itera- 
tion logic. 

"Ma'Cheirta tic ally , ■ divergence implies that the real world integrat- 
ing conditions acted in such a way that the guidance algorithm was 
unable to rectify the motion. Physically, divergence suggests that 
something in the dynamics or s/c system has been mismodeled, or under- 
designed, and that it has interactions with the other systems to cause 
wide, usually nonlinear, deviations from the reference mission. From 
the G&N point of view, it is these missions that are often of the 
greatest interest. They identify potential problems either in the 
baseline configuration or with the navigation and operational proce- 
dures. In many instances, divergence in the guidance can be traced 
to some well-understood phenomena, e.g., controllability, trajectory 
non-linearities, suboptimal schedule of guidance events, etc., but 
a clear identification and resolution of these problems in terms of 
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changes to the system design and/or operational procedures may not 
be so straightforward. 

The basic computational steps taken in a single iteration are 
as follows: First, an estimate of the actual state vector and the 

corresponding estimated integrating conditions are obtained from 
the simulated orbit determination logic. The estimated state is 
integrated to generate the estimated trajectory between the guid- 
ance point and the target. At the targeted stopping conditions, 
an estimated target error is computed by 

as * b. -h 

where T_ and T are the target variables on the estimated and refer 
ence trajectories, respectively * Next, a sensitivity matrix, S, of 
target variations with respect to control variations is computed 
about the estimated trajectory according to the matrix relation, 

S - q e 

where is the state to target transformation defined in Eq* 7-3 

and where 0 is that partition in the augmented state transition 
matrix which maps control changes into state variations (Eq. 7-7*5) , 
The matrix $ has the format. 


s 



T n> 


U ) 

m 


and has dimensions (nxm) , where n is the number of targets and m 
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the number of controls. With S It is possible to relate control 
changes at the maneuver to target variable changes according to 


L T - S A u. ' (7-12) 

If the matrix S is square, i.e., the number of controls and 

\ \ 

targets are equal, then the solution to (7-12) requires only a 
single matrix inversion, namely, 

A = S " 1 AT. (7-13) 

However, in most practical cases, £ is a nonsquare matrix (m>n) 
and a generalized inversion must be used, i.e. 

Am. = S T [sS T ] AT . (7-14) 


Note that Eqs. (7-13) and (7-14) again assume the form of a linear 
guidance relation A u = P AT . 

Once Au has been determined, the updated controls are used to 
generate a new estimated trajectory to see if the target errors have 
decreased. This overall process is repeated until the target errors 
are made less than some specified tolerances, or until a maximum 
number of allowable iterations has occurred. In SIMSEP, the prin- 
cipal measure of target error is given by a so-called quadratic error 
function defined by 


Q = 


.th 


■ £ 


i = i 


f ATQ) 

L t toi/^ 


sd 


where /^T(j) is the J component of the target error vector and 
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T T0L <j) is the component of a vector of target error toler- 
ances. Convergence occurs in the nonlinear guidance whenever Q is 
made less than one. 

Guidance divergence is said to have occurred if the quadratic 
error function is greater than the backup convergence criterion 
(A0K) after iterating a maximum number of times (NMAX) . Divergence 
also occurs if the quadratic error function increases on three 
successive predicted corrections. As far as the Monte Carlo mission 
currently being executed, divergence is considered to be catastrophic, 
and the mission is ended. As a backup, weak convergence occurs if 
strong convergence fails to be satisfied but Q is less than A0K. 

In this way, the nonlinear guidance algorithm can be tolerant of 
"near misses" without bringing the mission to a halt. 

Another feature included in the nonlinear guidance logic is 
the ability to weight certain low thrust controls more than others. 
This permits the user an added flexibility whereby he can effect a 
normalization of disparities in the units associated with controls. 

In addition, the user has the option to arbitrarily weight some 
controls more heavily, based upon knowledge and experience gained 
during the trajectory targeting process. 

Algorithms based on the computational steps outlined above are 
implemented in NLGUID. Both delta-velocity maneuvers and low thrust 
corrections are handled by essentially the same logic with A u in 
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Eqs. (7-13) and (7-14) being a velocity update for impulsive guid- 
ance. 

7.4 Simulated Orbit Determination 

SIMSEP, in the strictest sense of the word, is not a complete 
"simulation" in that an explicit orbit determination process is 
not included in the computational algorithms. The problem of 
estimating a state vector is done by sampling a knowledge covariance 
in much the same way as it samples a control covariance or an ephem- 
eris error covariance. Simply stated, an augmented knowledge 
covariance is sampled to obtain an error in the estimated state 
vector, £, relative to the actual trajectory state, Therefore, 

the estimated state vector is given by 


X 


5L + e 


Likewise, parameters which have been augmented to the state and 
estimated during the orbit determination process are also computed. 

Typically, the knowledge error covariances which are read as 
input to SIMSEP have been computed in an equivalent G0DSEP run, 
they are equivalent in the sense that the same trajectory and 
sequence of guidance events are evaluated. With this procedure, 
there is an added advantage of permitting a direct comparison of 
guidance results computed in G0DSEP and SIMSEP and their dependen- 
cies on the same state estimation results. 

There are several reasons why an explicit orbit determination 
capability has not been included in SIMSEP. Primarily, the 
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estimation process is not as subject to trajectory nonlinearities 

l 

as is the guidance process. This is because the estimation errors 
are generally small and well within the envelope of linearity. In 
addition, this method of simulating orbit determination minimizes 
the computational complexity of the program, while at the same time 
representing a cost effective means of performing an effective state 
estimation. 

7.5 Thrust Process Noise 

Digressing briefly to discuss the actual trajectory again, 
there is one very important process related to the generation of a 
real trajectory that is either ignored or modeled by an effective 
process in the other modes* This is the time-correlated thrust 
noise* These independent stochastic processes corrupt the commanded 
thrust controls, i*e., thrust magnitude, cone and clock angles, as 
small time-correlated perturbations. Each stochastic parameter is 
modeled in SIMSEP as a Gauss-Markov sequence which is computed during 
the actual trajectory integration. At time point the vector 

U , M of stochastic parameters is given by 



A — k + w 

“k+1 


where U, . has been evaluated at 


over the interval At = t, , 1 -t. , 

k+1 k 


t^* t*. is assumed to remain constant 

with its effect being determined by 


the coefficient matrix, A 
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where 0*^, T^» T n are correlation times associated with 

each corresponding stochastic parameter. —^+1 i s 3 vector of white 
noise terms which have statistics dictated by the requirement that 
the process remain stationary; namely 



■2 At/Tj ^ 2 
(1 - e ) 


where is the variance associated with the j component of 

the U k | vector. During the integration, ^ is evaluated at 

the start, the half-interval, and the end of a normal integration 


step. 

7.6 Guidance Execution Errors 

Once that a guidance correction has been formulated, the 
execution of that correction must be performed to affect the actual 
trajectory. The commanded correction computed by the guidance is 
an idealized set of control changes which are invariably corrupted 
by execution errors. For a low thrust control change, the executions 
errors are actually built into the thrust process through the thrust 
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biases and the dynamic process noise. However, for an impulsive 
maneuver, an explicit set of logic to corrupt the commanded delta- 
velocity change must be implemented. 

In general there are three basic execution errors which are 
modeled for impulsive maneuvers: 1) pointing, 2) resolution, 

and 3) proportionality. 'Given the commanded delta-velocity vector, 

A V . in its heliocentric representation the in-and-out of the 
ecliptic plane angles are determined as, 

U. = tan” 1 ( AV (2) / AV (3)) 
c c 

£ = Sin ” 1 ( AV C (3) /lAV^ ). 

Specified pointing angle errors are sampled to formulate changes 
in the commanded angles (S«<. and ) to simulate orientation errors for 
the actual delta -velocity, A V^ . Likewise, errors specified as 
proportional to the maneuver magnitude, Sp , and as a minimum meas- 
urable resolution of a maneuver magnitude, S f , are also sampled and 
added to the commanded correction. The actual velocity change is given as 

I A V A f = | Av { + S P + S r 

for the actual velocity magnitude. The actual velocity vector is 
given as 


L\a) = 


cos (^ + Set ) cos 

($+ S?) 

& V A P> - 

1 &Vjl 

sin(<t + Sot ) cos 


AV 3 > - 


sin(^ + S^ ) 
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APPENDIX 1 . 

t 

9 1 Conic Equations For Position And Velocity In E lliptical 
And Hyperbolic Orbits 

Given: r , v , t and y 

— O — T> O 

Find: r_ and v_ at time t. 

From the initial conditions we can find the inverse semi-major 
axis 


1 _ 2 



a r p 
o 

The mean angular motion 


n 


■ftp 


and relationships for the eccentricity and eccentric anomaly for 
elliptical orbits (a>o) 


e 


sin E 

o 



e 


cos E 

o 



or for hyperbolic orbits (a<o) 




r * v 

“O — o 

# sinh 

H 

o 

ii 

ii 



r v 2 

• cosh 

H 

o 

o o 
p 
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To find the position and velocity along the conic orbit. Keplers 
equation must be solved for the change in eccentric anomaly. 

Newton's method is used to solve the equation iteratively. Let 
Newton's method be given in the form 

t( V 

X k f 1 " x k f ' (x fc ) . 

Then for elliptical orbits 
x - E - E 

o 

f(x) = X + e-sin E (1 - cos x) - e-cos E sin x - n (t-t ) 

o o 

and for hyperbolic orbits 

x = exp (H - H ) 

° x 

f (x) = h e*exp(H o )-x + Ije-expC-lOx+l - ln(x+l) - n(t-t Q ) . 
The position and velocity at time t in elliptical orbits is given by 

r -{i - J o [ 1 - cos(E-E o >] } I, + i { sin(E - E o > 

- e* (sin E - sin E q )1 

- - Sin (E-E o ) • ^ i{l - f-[l - cos (E-E o )]l -V, 

o 

and for a hyperbolic orbit 


r-{l-i[cosh W-H 0 )-l]]r,+ « 

^e-(sinh) H-sinh H q ) - sinh (H-H)}- -v^ 

v = sinh (H-H ) r +{l - j[cosh (H-H ) -l] \ 

jl r r o ~~o L r u ° J 
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APPENDIX 2 

9.2 A Generalized 4th Order Runge-Kutta Algorithm With 
Runge 1 s Coefficients For A Matrix System Of First Order Diff erential 

Equations ♦ 

The 4th Order Runge-Kutta formula for numerically integrating 
first order Differential Equations of the form 

y* = f (x,y) 
is 

+ * * 2 + **, + v 


where h is the stepsize, x is the independent variable, y is the 
dependent variable and 

h ■ £,< vV 

V ^k*^l 

£ 2 - £ ' + r • y k + -r > 


\ 


^k’^2 

f 3 = f' (x k + - > 


f A = f’ (x k + h k , y k + V f 3> 

To generalize these equations for an mxn Matrix system of 
first order differential equations, write equation (1) as 


'k+l' 


(i, j) - y (i • j) + 5 s - [q (i.j) + 2-f 2 (i.J) + 2-fjU.j) + 
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and this can be written as 


\ 

Y... - Y. + ~ ( F 

k+1 k 6 1 


+ 2- F 2 i 2* F 3 + F 4 ^ 


where 


F, 


F' < V Y fc ) 


F 2 - F* (Xfc + 2 * Y k + 2 ‘ F l* 


F 3 =F-( V ),Y k 4.F 2 ) 


F 4 = F*(x k + , \ + \ ’.V 
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APPENDIX 3 

I 

9.3 Newton's 3rd Order Divided Difference Interpolation 
Polynomial . 

Given: (x^ , (x 2 , y 2 ) , (x^ y 3 > and (x 4> y^) 


Find: A third Order Polynomial that fits the given points 

Construct the following table 



where 


[ x i+i> *il = 


x i+l~ x i 


r l t x..+2 ,x i+il 

l x j+2’ X j+1» X jl x j+2 - x. 


, i = 1,2,3 


, j - 1,2 


l x k+3 * X k+2 , X k+l 


[ x k+3 ,X k+2 ,X k+ll t x k+2’ X k l ,X kl 


*k+3 " X k 


, k = 1 
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Newton's 3rd Order Divided Difference Interpolation Polynomial 
has the following form, 

y(x) = y^+(x-x) ♦ £x 2 ,x"] + (x-x^ • (x-x 2 ) • [ x 3 ,x 2 .x^ + 

(X-Xj^) * (X-X 2 ) • (X-X^) • [ .X^ ,X 2 jX^ 

while a 3rd Order Polynomial can be written 


3 2 

y (x) = ax + bx + cx + d 

To find a,b,c and d, the coefficients of the x terms in the Newton's 
Divided Difference Polynomial can be equated to a,b,c and d, such 
that 

a = [x 4> x 3 , x 2 , x x ] 
b = - (x 3 , x 2 , x^-a +[x 3> x 2> x x ] 

c = (x 3 , x 2 + x 3> x 2 ,x 1 )-a - (x 2 +x i )- [x 3 ,x 2 ,x 1 '] + [ x 2 » x il 

d - -(x 1 ,x 2 ,x 3 ) -a + x 1 x 2 [ x 3 ,x 2 ,x 1 ] - x x [x^x^+y 


Now a third order polynominal can be fitted to the four points and 
y can be determined from a given x or a maximum or minimum can 
be found from the following values of x, 


x = 


-b +s/b^ - 3 ^ 
3a 1 


is minimum 


and 


_ -b - - 

X 3a 


is maximum 
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APPENDIX 4 

9.4 Analytic Expressions for Terms in the Fa Matrix 

In Section 4.5 it was shown that the augmented state transition matrix, 
<p , is computed by integrating the matrix differential equation, 



In order to efficiently integrate this expression, it is necessary to have 
analytic representations for the individual elements of F^. 

has been identified in equation 4-6 as a matrix of first order 
partial derivatives obtained by expanding the equations of motion. In 
concise symbolism, may be written as, 

f a = .3 1&- 
dZA 

where f^ was defined in equation 4-5 and is the augmented dynamic state. 
For an analysis where covariances are propagated by the state transition 
matrix, i.e., the STM mode, the (maximum) component vectors of xa are 
defined as: 

S/C position vector 
S/C velocity vector 
thrust control vector 
harmonic coefficient 
Earth gravitational constant 
solar gravitational constant 


• r ' 


V 


u 

- 

J2 




CO 
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and the corresponding F^ matrix, in partitioned format, is 


0 I 


33 


L f 33 0 J66 

L 8 33J 63 

L k 3 1-* 6 1 

L d 3iJei 

L m 3li 

t° °J)6 

M 33 

*4 

Hi 

W 3 1 


Mb 

0 

0 

0 


[o] u 

0 

0 

0 

m 16 

Ml, 

0 

0 

0 


61 


(12x12) 

where the subscripts refer to the dimensions of the appropriate partitions, 


When the covariances are propagated by integrating the covariance differen- 
tial equation (see Section 4.6) in the PD0T mode, the augmented dynamic state 
vector is defined as 


f*. — 

r 


— — 
s/c position vector 

V 

= 

s/c velocity vector 

u 


constant thrust controls 

w 


time -varying thrust parameters (6x1) 


and F^ in this case is given as 


0 I 33 


'o’ 


"o' 


, f 33° 

66 

g 

63 

n 

66 

[0 0] 

36 

M 

33 

[«] 

36 

[° °] 

66 

[°] 

63 

[ h ] 

66 


(15x15) 
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In either STM or PDOT mode, the fully augmented state, as used by 
GODSEP, may also include measurement parameters. However, the terms 
appearing in F^ corresponding to measurement parameters are zero because 
they do not affect the dynamic process. 


Specific matrices in are defined as follows: 
f = partials of the s/c acceleration vector w.r.t. position 
components = d&.ldL j 


’33 


31 


31 


m 


31 


partials of the s/c acceleration vector w.r.t. the thrust 
controls - dji/dH , 

partials of S/C acceleration vector w.r.t. J2 in the gravitational 
potential function = 3 a/ 3J2, 

partials of the s/c acceleration vector w.r.t. the gravitational 
constant of the Earth = Ba/dfi^ > 

partials of the s/c acceleration vector w.r.t. the solar 
gravitational constant = da Jd[l s , 


n^g = partials of the s/c acceleration vector w.r.t. the time-varying 
thrust parameters = da/d w , and 

hgg = partials of the time derivative of the time-varying thrust parameters 
w.r.t. the time -varying parameters = q w/gw . 


As noted from its definition, the elements of the f^ matrix are evaluated 
by differentiating components of the S/C acceleration vector, £, with 
respect to S/C coordinates, i.e. 

f . .La 

33 2 L 



122 


where a is the sum of the n-body gravitational acceleration, g. , 
and the acceleration due to planetary oblateness, —J2‘ P ar ^^- a ^ s 

of with respect to r. are the components of the so-called gravity 
gradient matrix and are well-known, e.g., see Battin (Reference 4). 

In terms of S/C position vectors relative to gravitating bodies, r^ , 
the gravity gradient matrix is 


- I [3 1 , r. T - r . 2 I 33 ] fUr. 5 , 


with the summation being performed for all bodies. In this equation, 
yll . , is the gravitational constant of the i^ body, and I^ is a 
(3x3) identity matrix. 

The second factor comprising i s obtained by differentiat- 

ing the perturbing acceleration vector, a _ j2 , due to the nonspherical 
mass distribution in the primary. Differentiating each acceleration 
component with respect to x, y, and z generates a 3x3 matrix which is 
directly added to the 3^/9r matrix, evaluated above. The resulting 
matrix is given as 


3% 2 1 

^ — ^equatorial 


where 


" X 2 £ -D 

x y S' 

xzT 

xy$* 

y 2 S - u 

yzT 

xzT 

yzT 

z 2 T-3U 

15 J 2 M ^ 

11 ( 1 - 

N 

NJ 

T 7 

r 

V 

r 
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T 


15 

2 


J2 



( 



and 




This matrix must be rotated into an ecliptic representation for use 
in MAPSEP. This is accomplished by pre- and post-multiplication by 
the standard rotation matrix, E, which depends on the obliquity of 
the ecliptic. 


[ . E f 1%2 ~) e T 

*■ 3 ~ L 3E J equatorial 

The (3x1) k matrix is the matrix of partials which appear in the 
variational equations and reflects how variations in the magnitude of 
J 2 affect the S/C acceleration a. This matrix is given by 

9 I 

9 J 2 


where a^, a^, and a ^ are evaluated acceleration components given relative 
to the ecliptic coordinate system. 


0 

0 

0 


V J 2 

S y /J 2 

V J 2 



The 3x3 g matrix is given by 


g 


a y 



where A is the transformation matrix from spacecraft cartesian to inertial 

coordinates (Section 4.1) and transforms thrust controls to spacecraft 

coordina tes . 


0V' - 

a’ 

X 





— 

du 

a 1 

z 


-a 1 

sin (yaw) cos 

(pitch) 


a 1 

y 

0 


a 1 

cos (yaw) 



a 1 
z 

-a 1 

X 


a 1 

sin (yaw) sin 

(pitch) 

for the pitch/yaw system, with 

a T = 

thrust 

acceleration in spacecraft 

coordinates, and 






di' = 
9u 

r 

a * -a' 

x. y . 

-a' 

sinS 

cos Y 




a * a 1 

y x 

-a / 

sin S 

sinY 




a' 0 

z 

a 

cos S 





for the orbit plane system. 
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The 3x1 d vector is 



The 3x1 m vector is 



r 


r 


3 


The 3x6 n matrix is 



The 6x6 h matrix is 




where 


Te 


are process noise correlation times. 



Page 123-D has been deleted. 



Pages 124 through 128 have been deleted. 
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APPENDIX 5 


9.5 Tug Multiple-Impulse Orbit Transfer 

When the initial condition controls (r , r , i, f2 , 

P a 

w , f) are applied to TOPSEP, the reference parking orbit 
will change from iteration to iteration. In fact, it is pos- 
sible that the parking orbit characteristic of the last iteration 
may not be attainable directly from an Earth based launch within 

realistic launch constraints. Therefore, it becomes necessary to 
consider the interface between the SEP S/C and the launch vehicle 
in order to predict the "cost" of achieving the reference parking 
orbit. The cost may be estimated indirectly by determining the 
launch vehicle fuel budget to transfer the SEP S/C from some nominal 
inner parking orbit to the other reference parking orbit* If the 
inclination of the inner parking orbit is realistically constrained, 
an orbit plane change may be necessary to complete the transfer. As 
the angle of the plane change increases the estimated cost of the 
orbit transfer will increase dramatically. The estimated cost may 
then be used to distinguish between acceptable and unacceptable 
outer parking orbits while simultaneously sizing the fuel expenditure 
for presently conceived or operational launch vehicles (or intermediate 
stages) . 

The launch vehicle simplistically modeled in TOPSEP is the 
expendable space tug*. Initially, the tug is in a circular inner 


Although the launch vehicle is referred to as the space tug, in the 
remaining discussion, the specific vehicle is inconsequential to the 
sizing of the fuel budget (i.e., the vehicle may be a trans-stage. 
Burner II, Centaur, etc.) 
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inner parking orbit whose equatorial inclination is constrained 
due to bounds on the booster launch azimuth. The tug then 
performs the transfer to the outer parking orbit, which is the 
orbit specified by the initial state in the $TRAJ namelist. 

The tug's path from the inner parking orbit to the outer parking 
orbit will be a coplanar Holmann transfer if the required 
equatorial inclination does not violate the launch azimuth 
constraints. Otherwise, the tug follows a modified Hohmann 
transfer (i.e., a regular Hohmann transfer in the inner orbit 
plane followed by a plane change and circularization at the 
line of intersection with the outer orbit plane). The 
equatorial inclination of this transfer orbit is either the 
maximum or minimum inclination bound such that the required 
plane change is a minimum. For the coplanar transfer the 
ascending node of the inner orbit is fixed; thus, the launch 
azimuth may be computed explicitly and tested for a constraint 
violation. For the plane change transfer the launch azimuth 
is fixed, and the inner parking orbit is uniquely determined 
when the minimum plane change condition is enforced. 

Figure 9-2 illustrates the selection of the inner orbit 
normal vector which minimizes the plane change for the transfer 


to the outer parking orbit. 
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z 


ec 


normal vector of 
inner orbit mini 



Figure 9-2 Selection of the Normal Vector for the Inner 
Orbit Which Minimizes the Plane Change 


The larger cone in the figure represents the locus of normal vectors for 
all possible inner orbits having the maximum equatorial inclination 
allowed by the launch azimuth constraints (AZMIN and AZMAX in $T0PSEP). 
The smaller cone is the locus of normal vectors for orbits having the 
minimum equatorial inclination, which in general will be equal to the 
latitude of the Kennedy Space Center (28.608 deg). If the normal vector 
of the outer parking orbit falls between the two cones, the parking 
orbits are assumed coplanar. If the normal vector falls outside of this 
region, the inner parking orbit is characterized as follows: 
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(1) the inner orbit is inclined to the outer orbit by the 
angle where <p is the minimum angle between 

the normal vector of the outer orbit and the nearest 
cone, 

(2) the normal vector of the inner orbit becomes the 
projection of the outer orbit normal vector on the 
nearest cone. 

For both the Hohmann and modified Hohmann transfers two 

distinguishable velocity increments or Av 1 s will occur. Figure 9-3 

illustrates the relative positions where these maneuvers are executed 

for the general case (r specified by RP1 in $TOPSEP) . 

a 


z 

ec 
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The Hohmann transfer orbit is oriented so that the line of apsides 
corresponds to the intersection of the two parking orbits planes. 

The apoapsis radius of the transfer orbit is chosen to be the 
larger of the two radii to the outer parking orbit along the inter- 
section. The first impulse Av occurs at periapsis; the second 

— 3 

impulse Av occurs at apoapsis and provides a velocity increment 
lj 

to shape the orbit and to change orbit planes if necessary. If both 
parking orbits are coplanar the line of intersection becomes ambigious. 
For this case, it is assumed that the outer orbit and the transfer 
orbit share the same line of apsides. Thus r^ will correspond to 
the apoapsis vector of the outer parking orbit. 

The magnitudes of the velocity increments Av g and Av fa for 
the general case are computed as follows. 



where a is the semi-major axis of the outer parking orbit, 
o 

Since the first impulse does not initiate a plane change and since the 
velocity vectors before and after the impulse are aligned 

Av = (v ) + - (v )" 
a a a 
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The second impulse may perform a plane change through the angle <f> . 
Thus, by applying the law of cosines 

V \ 

Av b = ( (v b +)2 + (v b’ )2 " 2(v b +) (v b _) cos<p ) 

If <)> = 0, is the apoapsis radius of the outer parking orbit (r gQ ) 
and the above equation reduces to 

Av b ■ <V + - <V" 



where r is the periapsis radius of the outer parking orbit. 


Once Av and Av, have been calculated the fuel 
a b 

budget for the space tug may be computed. The tug's dry weight (W^ 
the propellent weight and tug's specific impulse (I 

are specified in the input namelist (TUGWT, TGFUEL, and TUGISP in 
$TOPSEP). The fuel required for the first maneuver (W^ is then 


) 

) 


J 


^fuePa 



W 


tot 



135 


where 


. = V/ + (W r J + w 

tot tug fuel max sep 


W p Is the weight of the SEP vehicle and g is the gravitational 
constant* It follows that the fuel for the second maneuver is 



The required fuel budget to perform the orbit transfers 
is then 


(W fueP 


= (W fuel>a + 


^fuelA 


tot 



Page 136 has been deleted 
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9,6 Control Weighting Schemes for TOPSEP 

Various weighting schemes are provided to allow the MAPSEP 
user flexibility in scaling controls in the TOPSEP mode (Chapter 5)* 
The purpose of these schemes is to alter the contours of constant 
cost and constant target error in the control space so that the 
projected gradient algorithm may converge more readily. Convergence 
problems occur most often when elements of the control vector differ 
in units. For example, there may be considerable difficulty in 
finding a converged solution when thrusting angles (cone or clock) 
are selected as controls in addition to thrust phase times. Whereas 
the internal units for the angles and times are radians and seconds 
respectively, the corresponding elements of the sensitivity matrix 
may vary by several orders of magnitude. That is, the sensitivity 
of the targets to a change of one radian in thrusting angle is many 
orders of magnitude greater than the sensitivity of the targets to 
a change of one second in thrust phase duration. In the example 
just described one would find that the PGM algorithm would compute 
a control correction which would try to eliminate the target errors 
by large changes in the thrusting angles and very small changes in 
thrust duration. Unfortunately, large control changes are often 
invalid in the very nonlinear control space associated with low 
thrust trajectories. To alleviate this problem, the normalized 
control weighting scheme has been devised. The diagonal elements 
of the control weighting matrix are defined as 
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Application of this weighting matrix to the controls determines a 
sensitivity matrix whose elements are roughly the same order of 
magnitude. Thus, the removal of the target error is spread evenly 
among the selected controls rather than among only a few* Hopefully 
all the control changes will be small enough to be valid in the 
nonlinear control space. 

Another type of convergence problem may occur. Sometimes 
elements of the sensitivity matrix vary by orders of magnitude 
even though the controls are all of the same units. For example, 
target parameters are much more sensitive to changes in thrusting 
angles early in the trajectory than they are to changes in these 
angles late in a trajectory. If the controls are not scaled, the 
PGM algorithm computes a control correction where changes in 
thrusting angles during early phases are unacceptably large and 
changes to the angles during later phases are undetectable. The 
following weighting schemes have been devised to spread the removal 
of target error more evenly among the selected controls, 
a) Sensitivity weighting 



b) 


Combined sensitivity, target error, and control weighting 



N 



i=l 


* 


u , 
J 


e . 
1 



Gradient weighting 



Averaged gradient and control weighting 



where G is defined in c 
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APPENDIX 7 


9.7 Integrated State Transition Matrices for Computing the 
Targeting Sensitivity Matrix 1 

Within the three basic modes of MAPSEP, trajectory guidance 
and/or retargeting represent one of the primary Computational 
problems worked in the program. Whereas the logic controlling 
these calculations is, in general, straight forward and easily 
understood, the actual execution remains as one of the more costly 
computational operations to be performed. This is especially true 
in TOPSEP and SIMSEP where targeting over long trajectory arcs is 
done repeatly. In order to minimize computational expenses, an 
algorithm which uses state transition matrices integrated with 
the trajectory has been implemented and is used exclusively in 
GODSEP and SIMSEP for computing the targeting sensitivity matrix. 
In TOPSEP, the user has the option to use either this or the more 
expensive, but equivalent, numerical differencing algorithm. 

The targeting sensitivity matrix, S, is a matrix of linear 
partials relating variations in the control variables, AH? to 
variations in the targets, A T , according to 

4 T = s A u • 

Looking at the targeting sensitivity matrix in more detail, it is 
seen to be of the form, 

^ (T^ , T^ T^) 

d (u 1 , ^2 > ***•» u m ) 


where the T's are selected target variables and the u f s are controls 
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Typical target variables include X^, Y^, Z^, r > v , etc, which are 
all evaluated about the final trajectory state. Typical low thrust 
control variables are the thrust phase stop time, thruster throttling, 
and thrust direction coefficients in the same or’ different thrust 
control phases. 

To provide a conceptual understanding of how S is evaluated 
from trajectory information generated in the augmented state transi- 
tion matrix, it is convenient to consider the trajectory segment 
depicted below where time points k. 



k + 1, ,,,, and f are shown bounding thrust control phases n, n + 1, 
etc. In the figure, time point f denotes the target time and represents 
the trajectory stopping condition. Considering a specific thrust con- 
trol phase, say n, it is recalled from Section 4-5 that the augmented 
state transition matrix generates partials of state variations of time 
k + 1 with respect to state changes at k and control variable changes 
interior to thrust phase n. In particular, these partials are contained 
in the | and 0 U partitions of the augmented matrix and may be*sym- 
bolically written as, 



142 


^(x, y, z, v x> v v> v z ) k + l ' 

<} (X, y , z, v x> X y , v z ) k 

y, z, v x , v y , v z ) k + x . 

^<V u 2 , u 3 , u 4 ) 

The u's correspond to the phase stop time, throttling, pitch (or 
in-plane) angle, and yaw (or out-of-plane) angle for the n phase. 
Other thrust policy coefficients are excluded from this treatment 
since their partials must be obtained by numerical differencing » 

If u^ (for example) in thrust control phase n is specified as 
an active control for the targeting event being considered, then the 
action of u^ on the state vector at f is computed by pre -mu It ip lying 
the appropriate column from ©^(n) for phase n by all intervening 
$ f s , i.e . 


and 


i - 


e = 


^ u 2 (n) 


1 


f, k+3 


1 


k+3, k+2 *k+2, k+1 


l 


^ . 

<iu 2 (n) 


Hence, an augmented say 9^, can be constructed by storing and 

pre-multiplying appropriate columns from the ® u ' s as they are computed 
during trajectory integration. The resulting gives the partials of 

the final state with respect to the active controls occurring in the 
various thrust phases between k and f and may be written as 


( x > y s 
^(u x , t 


v , 

X 


Vf 

U m ) 


u 


9 
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The final component necessary to compute the requisite S matrix 

is the evaluation of target variable partlals with respect to the 

final state, L- This is most expediently done by a numerical dif- 
— r 

ferencing algorithm to generate the differential point transformation 
matrix, '•'j . The matrix can be written as 

= i (T i» T 2 a •"» V . 

\ 4 (*» y, z, v V Vf 


Now S is seen to be 


V- 


In TOPSEP where initial state conditions are also permitted as 
control variables, it is necessary to extend the above procedure but 
not the computational method. For this special case, only the chained 

I ’s are needed to compute the partials of final state with respect 

* | 

to changes in the state at k. Clearly, selected columns from 
which is defined by 

If , k If, k+3 Ik+3 , k+2 Ik+2 , k+1 Ik+1 , 


may be augmented to the previously defined 9^ to give the requisite 
targeting sensitivity matrix to be used by TOPSEP, 
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APPENDIX 8 


9.8 The Shadow Model 

Passage into the shadow of the Earth and the occultation time 
are especially important considerations for Earth orbital SEP missions. 

The occultation of the sun, which has the effect on the S/C of an 
imposed coast period, restricts the mission thrusting time. Additionally, 
thruster efficiency is reduced due to cyclical cooling and warm-up 
periods for the thruster subsystem. In this appendix, the MAPSEP 
shadow model is developed. A general technique for determining shadow 
entrance and exit times is reviewed (Reference 10) and the thruster 
warm-up model is discussed. 6 

The analysis of eclipse times is conveniently simplified by the 
following assumptions: 

1. No flattening of the Earth; 

2. No shift of the Earth in its orbit during occultation 
of the sun; 

3. No umbra and penumbra effects (cylindrical shadow); 

4. No orbital changes due to thrust between the time the 
shadow computations are completed and the predicted time 
of shadow entrance. 

The error which is introducted by the latter assumption is minimized by 
forcing the shadow computations to occur immediately preceding the 
estimated time of shadow entrance. This is in fact the process enforced 
within the program trajectory propagator. For additional information 
concerning implementation of the shadow model, refer to subroutine PATH 
in the Program Manual . 
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The preceding assumptions fix the orbital geometry in the ecliptic 
reference frame so that it is a simple matter to determine whether the 
orbit intersects the cylindrical shadow. Once it has been established 
that an intersection exists, a quartic equation in the cosine of the 
entrance and exit true anomalies may be formulated. The corresponding 
shadow entrance and exit times may be computed using Kepler's equation. 

The elliptical orbit is completely defined as to size and orienta- 
tion by the classical orbital elements, 

a, e, i, fl , a, f 

or equivalently 


a, 


p, Q> 


f 


where P and Q are the unit position vector and unit velocity vector 
respectively at the point of periapsis. In terms of i, ft , co the 
unit vectors may be defined as 


A 

P = 


and 


3 = 


cos co o cos ft -sin co • sin ft ° cos i 

cos co • sin ft -hs in to * cos ft • cos i 

sin co * sin ft 


-sin co * cos ft -cos co * sin ft < cos i 
-sin co • sin ft +cosw * cos ft ° cos i 
cos co • sin i 


9 , 8.1 


9 . 8.2 
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If the orbit is circular, P is arbitrarily defined as the unitary node 
vector in the ecliptic ref erence .frame . For the special case when the 

A 

orbit is circular and in the ecliptic plane, P is defined as a unit 
vector in the X-direction. 

It is not necessary that the orbit be explicitly defined by the 
complete set of orbital elements before testing the occultation^ however* 
The following coarse test can be made to determine if the orbit which is 

A A 

partially defined by the elements a, e, and N intersects the shadow* N 

A A A 

is simply the unitary orbit normal vector (i.e., N - P x Q) . The objec- 
tive of this test is to eliminate further shadow computations if the 
orientation of the orbital plane is such that it is impossible for any 
ellipse of size, a, and of shape, e, to intersect the shadow. 
Certainly, if the point of periapsis, does not fall within the 

A 

shadow given any orientation of P in the orbit plane, then occultation - 
of the sun cannot occur. The geometric relationships are obtained from 
Figure 9.8*1. 



Figure 9.8.1. 


Orbit Geometry 
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Consider the angle 
sin a = 

and the critical angle 
sin a 


0 (See Figure 9.8.2) which is defined 

A A 


N • S 


O which is defined 
cr 


= R /r 
cr e p 


9.8.3 


9.8.4 


If 

| sin a | > |sina cr | , 9,8,5 

then no intersection can possibly exist between the orbit and the shadow* 
However, if the relationship in equation 9*8.5 is not true, a definitive 
conclusion cannot be reached and further tests must be made* 

If P and Q are specified in addition to a and e, conclusive results 
can be obtained. First, the projection of S on the orbital plane must 
be found. Let T be the transformation matrix defined T = [ P | Q ] • 
Then the projection of S on the orbital plane is 


s 


A. 

S 


9.8.6 


and the true anomaly of s. in the orbit is 


$ = tan 1 (s 2 / s 1 ) 


The position vector magnitude at 0 in the orbit is 


9.8.7 
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r = p/(l + e cos Q ) 9.8.8 

s 

where p is the semi-latus rectum. Again consider the angle 0 which 
is defined in equation 9.8*3 and the critical angle which is defined 


sin 


a 


cr 


R fr 

e s 


9.8.9 


If the relationship from equation 9.8.5 is true, no Intersection can 
possibly exist between the orbit and the shadow. If the relationship 
is false, it is assumed that part of the orbit falls within the 
cy 1 i ndr ic a 1 shad ow . ^ 

The eclipse time can be estimated in terms of the angles displayed 
in Figure 9.8.2. This time can then be used to identify the set of 
osculating orbital elements from which the quartic equation is to be 
formulated . 



Figure 9.8.2. Shadow Geometry 

^“This test may be insufficient for highly elliptical orbits with the line 
of apsides nearly aligned with *S and 0=(X cr . If this situation occurs 
and equation 9.8*5 is false, an intersection may not exist. 
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The angles A0, a, and (J cr are related by the equation 

COS ( (X ) = COS ( O ) cos (A/3 ) 
cr 


so that 


A0 


-1 

COS 


(cos a cr / cos a 


) 


9.8.11 


The transit angle through the shadow is simply 2A/J . From the 
conservation of angular momentum equation 


Af h ( 1 + e cos f) 
r At “ -P 


Solving for At » 


At 


Af 


(1 + e cos f) 


2 



9.8.13 


where h and fJ are the angular momentum and planetary mass respectively. 
From equation 9.8.13, the shadow time may be estimated. That is^ 


t 

s 


2 A/2 

(1 + e cos /3 



9.8.14 
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Thus, the time of shadow entrance, t^, relative to periapsis crossing 
can be computed from the eccentric anomaly of s_ in the orbit and the 
estimated shadow time. That is, 

E = 2 tan ^ (tan B /2) / 1 " 9.8.15 

s yj 1 + e 


so 


t = E - e sin E - t /2. 9.8.16 

q s s s 

In order to predict the shadow entrance and exit times more 
accurately, a quartic equation in the cosine of the entrance and exit 
true anomalies may be solved. To minimize the errors due to assumption 
4., the osculating orbital elements are selected at time t^. The 
following vector relationship is obtained from Figure 9.8.1. 

9.8.16 


9.8.17 


d + R = r 


Upon entrance to or exit from the shadow 


A 

S * d = d 


or 

S • (r - R e ) = d 


9.8.18 
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since S is perpendicular to R g . Therefore, the angle between the anti 
sun vector and the radius vector to the SEP S/C at entrance or exit is 


2 

cos (// = (r^ - R 0 ) ^ / r 

■* 7 T \ ijj \ 77 * 

The angle is taken to be in the interval ~ ? S' ? 2 

since this is the only domain where a shadow exists. At all times, 
the angle ip can be obtained from the equation 


cos ip - (s«r)/r 

but 


r - T 

r cos f 


r cos f 


9.8.20 


9.8.21 


which illustrates the mapping of the radius vector to the S/C from the 
orbital to the ecliptic coordinate system. Thus, 

cos ip = (r cos f) S«P + (r sin f) S«Q 9,8.22 

r 

If the dot products in equation 9.8.22 are defined as 

A A- _ A A 

A - S»P and P - S 'Q 

where the coordinates of S, P, and Q are taken at time t q , there results 
the simplified equation 
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cos {p = A cos f + Psin £ 9,8.23 

When the squares of equations 9.8.19 and 9.8.23 are equated and the 
relation r = p / (1 + e cos f) is employed, the shadow function is 
produced: 

- (1 + e cos f) 2 R^ 2 + p 2 ( A cos f + Psin f) -p 9.8.24 

where = 0 is the condition for entrance or exit from the shadow. 

Note that only those solutions of equation 9.8.24 obtained when 

A cos f + psin f y 0 9.8.25 

are of any physical meaning. If equation 9.8.25 is written in the form 

= A cos^ f + A cos 3 f + A n cos 2 f + cos f + A 

O 1 £ -5 3 

• where 

A o = (Re/p) 4 e 4 - 2 (Re/p) 2 ( p 2 - A 2 ) e 2 + (p 2 + A 2 ) 2 

A^ = 4 (Re/p) 4 e"^ - 4 (Re/p) 2 (p 2 - A ) e 

A 2 = 6 (Re/p) 4 e 2 - 2 (Re/p) 2 ( p 2 - A 2 ) - 2 (Re/p) 2 (1 - p 2 )e 2 

+ 2 (p 2 - A 2 ) (l - p 2 ) - 4 p 2 A 2 


9.8.26 
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A 3 = 4(Re/p) 4 e - 4(Re/p) 2 (1 - p 2 ) e 

A 4 = (Re/p) 4 - 2 (Re/p) 2 (1 - p 2 ) + (1 - pV 

2 

then equation 9.8.24 is solvable in closed form by quadratic radicals. 
Rejection of spurious roots is accomplished by application of equations 
9.8.24 and 9.8.25. Once the entrance and exit true anomalies are formed, 
the eccentric anomalies may be computed as in equation 9.8.15. If ^ nter 
and E are the entrance and exit eccentric anomalies respectively, then 


and 


t _ = E ■ e sin E 

enter nter nter 


'exit - E xit ' e sin E xit 


and the total time in the shadow is then 


t s fc exit t enter 


9.8*28 


Whenever the SEP S/C leaves the shadow, it must initiate a warm-up 
period before it can operate at full efficiency. For any given shadow 
time, the engine-restart delay is obtained from the function 


0, t s < c 

2 

a + a. t + a t ,t > C 
o Is 2 s * s ' 


9,8.29 


A closed-form solution to quartic equations may be found in Escobal's 
Methods of Orbit Determination . "Appendix III, Pages 430-434. 
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D is considered to be zero if the time in the shadow is less than C 
Figure 9.8.3 illustrates the warm-up delay model when a Q = 9, 

= .464, a 2 = 0 and C = 3 (e.g., the default parameters in the 

$TRAJ namelist) . 



Shadow Time (min) 


Figure 9*8,3 Warm-Up Time 




